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Perturbation  methods  were  used  to  study  the  propagating  fields  inside  a 
weakly-guiding  circular  fiber  whose  refractive  index  profile  is  slightly  deviated  from 
the  optimal  truncated-parabolic  shape.  It  was  observed  that  the  performance  of  the 
fiber  is  reduced  dramatically.  Next,  finite  element  methods  are  considered  to  solve 
Sturm-Liouville  and  Helmholtz  equations  defined  in  unbounded  domains  so  that  the 
non-physical  solutions  can  be  eliminated  from  the  computations.  It  is  known  that 
a circular  fiber  with  a smaller  core  yields  more  evanescent  energy.  From  numerical 
results,  it  was  observed  that  a step-index  square  fiber  whose  core  cross-sectional  area 
is  slightly  larger  than  that  of  a step-index  circular  fiber  still  has  more  evanescent 
energy  than  the  circular  fiber. 
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CHAPTER  1 
INTRODUCTION 


Optical  fibers  have  many  important  applications  in  telecommunications, 
evanescent  wave  sensoring  and  image  processing.  There  is  a vast  amount  of  research 
papers  and  books  in  electrical  engineering  and  applied  mathematics  devoted  to  the 
study  of  their  properties  and  performance.  In  this  work,  we  present  some  of  the 
approximation  methods  that  are  useful  in  estimating  the  propagating  electromagnetic 
fields  in  circular  and  non-circular  fibers. 

In  telecommunications,  circular  weakly-guiding  fibers  are  used  to  transmit 
information  of  several  hundred  Mbits/sec.  They  are  characterized  by  the  small  dif- 
ference in  the  core  and  cladding  refractive  indices.  In  this  case,  Maxwell’s  equations 
can  be  well  approximated  by  Helmholtz  equations  by  ignoring  the  polarization  effect 
(or  first-order  approximattion).  These  equations  can  be  further  reduced  to  Sturm- 
Liouille  equations  using  polar  coordinates.  This  simplication  is  possible  due  to  the 
circular  symmetry  of  the  fibers.  It  has  been  shown  that  multi-mode  fibers  with  a 
truncated  parabolic  refractive  index  profile  achieve  the  optimum  bandwidth-the  in- 
formation carrying  capacity  [5].  In  the  process  of  fiber  fabrication,  however,  it  is  not 
possible  to  obtain  the  exact  parabolic  shape  for  the  refractive  index  profile.  There- 
fore, it  is  of  importance  to  examine  the  effect  on  the  fiber  bandwidth  when  its  index 
profile  is  slightly  deviated  from  the  ideal  shape.  Perturbation  methods  are  often  used 
for  this  type  of  problem.  It  is  known  that  a small  perturbation  to  the  optimal  profile 
reduces  the  fiber  bandwidth  dramatically  [30,  15,  24,  25].  This  negative  effect  can  be 
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predicted  by  calculating  the  Fourier  transform  of  the  impulse  response  of  the  fiber 
which  is  defined  in  terms  of  propagation  constants  /?.  At  the  present,  many  fiber 
manufacturers  still  use  the  first  order  perturbation  method  developed  by  Olshansky 
[30].  In  this  method,  the  transverse  electric  field  for  the  optimal  profile  is  approx- 
imated by  the  Laguerre  solution  for  the  wave  guide  with  infinite  parabolic  profile 
and  the  study  of  field  perturbation  is  carried  out  in  terms  of  the  Laguerre  functions. 
Marcuse  et  al  later  presented  a different  approach  to  the  same  problem  using  WKB 
approximation  [24,  25).  However,  this  method  is  known  to  fail  at  turning  points. 
We  wilt  present  a better  approximation  using  first  and  second  order  perturbation 
methods  based  on  the  exact  analytical  solutions  (eigenvalues  and  eigenfunctions)  of 
the  optimal  profile.  These  exact  solutions  can  be  expressed  in  terms  of  Whittaker 
functions  of  the  first  kind  and  modified  Bessel  functions  [32,  34].  We  also  give  the 
estimates  for  the  coefficients  of  the  Taylor  expansion  of  the  perturbed  eigenvalues. 

For  arbitrary  index  profile  shape,  we  rely  on  finite  element  methods  to  find 
the  solutions  of  singular  Sturm-Liouville  equations.  The  singularity  arises  from  the 
unboundedness  of  the  considered  domain.  As  a result,  weighted  Sobolev  spaces  and 
modifed  finite  element  spaces  are  needed.  In  this  case,  the  usual  finite  dimensional 
approximating  spaces  are  replaced  by  infinite  dimensional  spaces  to  reduce  the  com- 
putational domain  and  to  eliminate  non-physical  solutions  called  spurious  modes. 
This  technique  is  discussed  in  [17]  where  an  unbounded  interval  is  reduced  a bounded 
one,  /,  on  which  the  usual  finite  element  functions  are  used  {e.g.,  piecewise  linear 
functions).  Outside  of  /,  all  functions  with  compact  support  in  1R\  / will  be  con- 
tained in  the  approximate  subspaces.  This  method  yields  the  same  convergence  rates 
for  eigenfunctions  and  eigenvalues  as  in  the  standard  bounded  interval  case,  i.e.,  the 
order  of  the  eigenvalue  error  is  the  square  of  the  order  of  the  interpolation  error  if  the 
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eigenfunctions  are  sufficiently  smooth.  To  use  this  method,  however,  one  has  to  have 
the  fundamental  solutions  at  each  end  of  the  bounded  interval  / as  boundary  func- 
tions. In  our  problem,  only  the  fundamental  solution  at  infinity  is  available.  Thus,  it 
is  natural  to  generalize  this  method  by  using  appropriate  approximating  subspaces. 
We  also  obtain  the  same  convergence  rates  for  eigenvalues  and  eigenfunctions. 

If  the  core  of  the  fiber  is  not  circular  but  rather  rectangular  or  of  any  other 
shape,  Helmholtz  or  Maxwell’s  equations  have  to  be  solved.  A numerical  method  de- 
veloped in  [16]  uses  the  circular  harmonic  series  representation  of  the  electromagnetic 
fields  in  rectangular  fibers,  f.e.,  the  series  of  Bessel  and  modified  Bessel  functions. 
Though  one  can  obtain  a solution  as  accurate  as  desired  by  increasing  the  number 
of  terms  in  the  series  expansion,  the  computation  is  complicated  and  costly.  A more 
attractive  alternative  method  is  given  in  [23].  In  this  method,  the  exterior  corner 
regions  of  a rectangular  waveguide  is  completely  ignored.  This  approximation  is  jus- 
tified only  when  the  mode  is  not  very  close  to  cut-off.  However,  in  fiber  optic  sensor 
application,  the  modes  near  cut-off  play  an  important  role  since  a large  fraction  of 
their  energy  leaks  into  the  exterior  region  which  can  be  used  for  chemical  detection. 
Thus,  finite  element  methods  again  become  the  obvious  choice  for  solving  the  problem 
numerically.  The  electromagnetic  fields  of  guided  modes  in  a fiber  satisfy  Helmholtz 
equations  in  the  whole  IR^.  Since  the  domain  is  unbounded,  it  is  custom  to  impose 
the  zero  Dirichlet  condition  u = 0 on  an  artificial  boundary  sufficiently  far  away 
from  the  fiber  core.  This  method  is  justified  by  the  exponentially  decaying  property 
of  guided  modes  outside  the  core.  The  portions  of  guided  modes  that  exists  in  the 
cladding  is  called  evanescent  waves.  They  have  applications  in  spectroscopy  and  fiber 
optic  evanescent  sensoring.  Unfortunately,  this  standard  method  is  known  to  have 
many  spurious  (non-physical)  solutions.  This  phenomenon  has  been  observed  and 
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studied  by,  for  example,  Rappaz  in  [31].  To  eliminate  these  spurious  solutions,  we 
borrow  a technique  known  as  Dirichlet-to-Neumann  maps  (DtN  maps)  developed  for 
exterior  and  scattering  problems.  This  technique  introduces  an  artificial  boundary 
and  imposes  on  it  a relationship  (map)  between  Dirichlet  and  Neumann  boundary 
conditions.  A good  reference  to  this  useful  method  is  [14].  The  Dirichlet-to-Neumann 
maps  can  be  applied  successfully  to  our  unbounded  problem.  It  also  turns  out  that 
the  error  estimates  are  easily  proved  using  the  positive  definite  property  of  DtN  maps 
and  a trace  theorem.  We  apply  this  method  to  calculate  and  compare  the  evanescent 
energy  of  a step-index  square  fiber  and  that  of  a step-index  circular  fiber.  This  energy 
existing  outside  of  the  core  of  a fiber  is  useful  in  optical  spectroscopy.  For  numerical 
experiments,  we  choose  a step-index  square  fiber  with  a slightly  larger  cross-section 
area  than  the  cross-section  area  of  a step-index  circular  fiber.  This  choice  is  justified 
by  the  fact  that  in  general  a fiber  with  a smaller  core  has  more  evanescent  energy.  In 
spite  of  having  the  larger  core,  the  square  fiber  is  observed  to  have  more  evanescent 
energy  existing  outside  of  the  core  than  the  circular  fiber. 


CHAPTER  2 

WEAKLY-GUIDING  APPROXIMATION 
2.1  Circular  Optical  Fibers 

An  optical  fiber  is  a three  dimensional  cylindrical  medium  consisting  of  two 
major  parts:  the  inner  core  made  of  dielectric  material  like  silica  or  composite  mate- 
rial of  different  refractive  indices  (graded  fibers)  and  the  outer  cladding  surrounding 
the  core,  usually  made  of  a lower  refractive  index  material  like  plastic  or  silica  non- 
doped.  The  cladding  is  considered  to  be  extended  to  infinity  in  the  radial  direction 
for  purposes  of  analysis.  This  assumption  is  justified  in  practice  where  the  cladding 
is  much  larger  than  the  core  in  diameter.  We  orient  the  fiber  so  that  its  longitudinal 
axis  is  parallel  to  the  z-axis  in  IR^.  We  also  assume  that  the  geometry  of  the  trans- 
verse cross-section  and  refractive  index  distribution  of  the  fiber  are  2:-translationally 
invariant,  Le., 

n{x,  y,  z)  = n{x,  y),  (x,  y)  G IR^ 


Figure  2.1:  Optical  fiber 
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Let  0 X IR  denote  the  core  region  where  0 is  the  two  dimensional  transverse 
cross  section.  Denote 


rico  = max  n(x,y), 
{x,y)en 

rici  = cladding  index. 


(2.1) 


In  order  for  a guidance  to  take  place  in  a fiber,  the  guiding  condition  must  be  satisfied, 
i.e., 

rico  > rich  (2.2) 


Definition  2.1.1  The  normalized  index  difference  A is  defined  as 


^ ^ 1 


ni 


(2.3) 


This  parameter  is  important  in  fiber  optics.  There  are  two  major  classes  of  commer- 
cial fibers:  multi-mode  and  single-mode.  A single-mode  fiber  allows  one  propagating 
mode  while  a multi-mode  fiber  has  more  than  one  propagating  mode.  A multi-mode 
fiber  usually  has  the  core  radius  between  25  and  50//m  and  the  normalized  index 
difference  A « 1%  . A single- mode  has  the  core  radius  between  2 and  5 //m  and 
the  normalized  index  difference  A « 0.2%.  The  standardized  measurement  for  the 
cladding  diameter  is  125  /xm. 

The  electromagnetic  fields  in  a fiber  are  governed  by  Maxwell’s  equations. 

2.2  Time  Harmonic  Maxwell’s  Equations  and  Transmission  Conditions 

The  general  Maxwell  equations  for  electromagnetic  fields 
(E(x,y,z,f),H(x,2/,z,f))  are 


dt 


V X H = 0, 


^ ^ X E — 0. 


(2.4) 
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where  /tq  is  the  free-space  magnetic  permeability,  e (x,y,z)  = Cq  n'^{x,y,z)  is  the 
electric  permittivity  of  the  medium  , Eq  is  the  free-space  electric  permittivity,  and 
n is  the  refractive  index  distribution  of  the  material.  In  our  treatment,  waves  are 
assumed  to  propagate  in  linear  media,  that  is,  e is  independent  of  field  strength. 

For  a time  harmonic  wave,  its  electromagnetic  fields  are  expressed  as 
E{x,  y,  z,  t)  = E(x,  y,  z)e~''^\ 

(2.5) 

H(a:,  y,  ^r,  t)  = IH(a;,  y,  z)e  ‘"h 

where  cj  is  the  angular  frequency.  Let  c be  the  speed  of  light  in  free  space,  then 
c = and  k = ^ ^ is  the  wave  number. 

Time  harmonic  Maxwell’s  equations  are  of  the  form 

V X E = io;/xoIH, 

V X IH  = — iweE, 

(2.6) 

V • (eE)  = 0, 

V • H = 0. 

V 


Transmission  conditions.  When  the  electromagnetic  field  crosses  a surface 
S which  is  a boundary  of  two  distinct  media  of  different  refractive  indices,  it  is 
necessary  to  impose  the  following  conditions  known  as  transmission  conditions.  Let 
n denote  the  normal  unit  vector  to  S pointing  from,  say,  medium  1 to  medium  2.  Let 
E6),H(‘),i  = 1,2,  be  the  electromagnetic  field  in  the  medium  i,i  = 1,2.  Then  on  S 
we  impose  the  following: 


E(i)  X n = E<2)  X n. 

(2.7) 

X n = 1H(2)  X ri. 

(2.8) 

?E(i)-n  = 

(2.9) 

lR(‘^-n  = lH<2)-ri, 

(2.10) 
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where  ui  and  U2  are  the  refractive  indices  of  medium  1 and  medium  2,  respectively. 
In  the  case  of  optical  fibers,  the  surface  S is  usually  a circular  cylinder  and  the  unit 
normal  vector  h = (nj:,ny,0). 

Definition  2.2.1  A guided  wave  is  a particular  solution  (IE,  H)  to  time  harmonic 
Maxwell’s  equations  (or  just  Maxwell’s  equations  ) such  that  IE  and  IH  satisfy  the 
transmission  conditions  (2.7-2.10),  of  the  form 

E = E{x,y)e~'^\ 

IH  = 

where  E G [L^(IR^)]^  and  H G [Z/^(IR^)]^,  i.e., 

f (|E;|2  + |/f|2)  dxdy<  oo, 

where  (d  > 0 is  the  propagation  constant  and  kn^  < fd  < knco  ■ 

The  definition  of  guided  waves  signifies  the  fact  that  these  waves  are  plane 
waves  propagating  without  any  distortion  in  the  fiber  axis  direction  (^-direction).  The 
vector-valued  two  dimensional  vector  fields  E{x,y)  = (Ex{x,y),  Ey{x,y),  E^{x,y)) 
and  H{x,y)  = {Hx{x,y),  Hy{x,y),  H^{x,y))  describe  the  the  distribution  of  electro- 
magnetic field  in  each  cross-section. 

2.3  Weakly-Guiding  Approximation 

Most  of  commercial  optical  fibers  used  in  telecommunications  are  of  weakly 
guiding  type,  i.e.,  the  normalized  index  difference  A is  very  small.  Consequently, 
Maxwell’s  equations  can  be  well  approximated  by  a simpler  set  of  2D  scalar  wave 
equations  known  as  Helmholtz  equations.  These  equations  can  be  further  simplified 
to  ID  equations  if  the  geometry  of  the  cross-section  is  a circular  disk.  In  this  section. 
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we  derive  the  Helmholtz  equations.  The  \D  equations  will  be  addressed  in  the  next 
chapter.  We  first  return  to  Maxwell’s  equations  (2.6).  From  the  first  two  equations 
we  can  eliminate,  for  example,  the  vector  field  IH  to  obtain  an  equation  for  IE  only, 
i.e., 

V X VIE  = tuVo  £ IE  = Ar^n^IE.  (2.11) 

Using  the  vector  identities 

V X V-JF  = - V(V  • JF), 

V X (gT)  = gV  ■ -Vg, 

where  .F"  is  a vector-valued  function  and  g a scalar  function,  we  have  in  each  region 
with  smooth  that 


V^E-l-kVlE  = V(V-E) 


= -V(Et  ■ Vt  ln(n2)),  (2.12) 

where  Et  :=  Ete~‘^*  ;=  (Ex,  Ey)e“‘^*  is  the  transverse  component  of  E and  V(  := 
(^,  ^)  is  the  transverse  component  of  the  differential  operator  V. 

Since  E :=  (Et,  Ez)e~‘^®,  we  can  also  write  (2.12)  as 

( = _(Vt  - i/?z)(Et  • Vt  In(n^)),  (2.13) 

where  z — (0,0,1)  is  the  unit  vector  in  ^-direction.  Then  the  equation  for  the 
transverse  component  of  E is 


(V?  + kV  - = -V,(E,  ■ ln(n")). 


(2.14) 
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The  term  V(ln(n^)  represents  the  polarization  effect  that  exists  when  is  a non- 
constant function.  It  is  convenient  to  express  n^(x,y)  as 

n^(x,  y)  = nl^{l  - 2Af{x,  y)),  f{x,  y))  > 0.  (2.15) 

We  then  expand  in  Taylor  series  in  terms  of  A and  obtain 


V,ln(n2)  = V*ln(l-2A/(x,y)) 

= V,(-2A/(x,y)-2AV"(x,j/)...) 

= -2AVtf{x,y)-2A^Vtf{x,y)-... 

If  A 4;  1,  the  first-order  approximation  of  (2.14)  is 

(V? -H  fcV  - ^'*)IEt  = 0.  (2.16) 


This  approximation  is  equivalent  to  ignoring  the  polarization  effect  and  is  known 
as  weakly-guiding  approximation.  The  equation  (2.16)  is  known  as  a Helmholtz 
equation. 

Definition  2.3.1  Assume  A 4C  1.  The  weakly-guiding  approximation  of  (time  har- 
monic) Maxwell’s  equations  in  terms  o/IEt  is  of  the  form 

(V2  + Pn^  - = 0, 


or  equivalently, 


(V^-f  A:V-/?2)£;t  = 0 


(the  factor  e is  omitted). 

In  a similar  fashion,  we  can  derive  the  weakly-guiding  approximation  for  IH.  From 
the  first  two  equations  in  (2.6)  we  have 

V X (-V  X H)  = 
e 


(2.17) 
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Using  the  vector  identity 

V • {9^)  =9'^  X T +Vg  • 

we  obtain 

-^V  X V X IH  - V(l/n^)  X e = k^H, 
thus,  the  transverse  component  of  IH  satisfies 

(V,"  + - /3^)mt  = Vt  InCn"*)  X IH.  (2.18) 

Hence,  the  weakly-guiding  approximation  yields 


(V?  + fcV  - /?2)IHt  = 0.  (2.19) 

Definition  2.3.2  Assume  A 1.  The  weakly-guiding  approximation  of  Maxwell’s 
equations  in  terms  of  IH(  is  of  the  form 

(V2  + ^2^2  _ ^ 


or  equivalently, 


(e'^”  is  omitted). 


( y2  + - p^)Ht  - 0 


The  longitudinal  components  and  of  E and  H respectively  can  also  be  shown 
to  satisfy  the  scalar  Helmholtz  equation  under  the  weakly  guiding  approximation. 
By  writing  out  the  equations  V x IE  = iu;/zoIH  and  V x IH  = — icuelE  we  get  a system 
of  six  equations: 


dE, 

dy 


-|-  iPEy 


-iPE,;  - 


dE, 

dx 


iu)p.Q  Hx , 

iujpoHy, 
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dEy  dE^ 


dx 

dH, 

dy 


dy 


+ i(3Hy  - —iujeEx, 


-i(3Hx  - 


dH. 


dx 

dHy  dH^ 


= —iujeEy, 

- —iueEz. 


dx  dy 

From  these  equations,  we  can  solve  for  E^,  Ey,  Hx,  Hy  in  terms  of  the  longitudinal 
components  E^,  and  their  partial  derivatives,  i.e.. 


Ex  = 


Ey  = 


Hx  = 


,.dEz^  dHz. 


k'^nP'  — dx  ' dy 
i ^JEz  dHz 

dy 
i 


k'^n?  — dx 


[p— u;£o« 


dx 

dE 


), 


dy 


), 


dy 


^^dHz  , 2 9Ez^ 

[P-^  +o;£on  yo-^)- 


Remark:  If  we  set 


E = (Ex,Ey,-iEz), 

H = {Hx,Hy,iHz), 

then  this  modification  will  allow  us  to  work  in  the  real  space  IR^  instead  of  C?. 
Thus,  it  is  enough  to  solve  for  the  longitudinal  components  Ez  and  Hz  of  E and  H 
respectively.  The  ^-component  of  V x V x IE  = k^n^IE  satisfies 

+ ln{^)  (2.20) 

= f/?£^<Vln(n^),  (2.21) 

and  the  ^-component  ofV  x (;^)V  x IH  = k^IH  satisfies 


V^Hz  + (fcV  - f3^)Hz  = -{VHz  - i/3Ht)Vln{n^). 


(2.22) 
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As  before  , if  we  ignore  the  polarization  effect,  we  obtain  the  scalar  Helmholtz  equa- 
tions for  E:s  and 


-f  {k'^n^  - (i^)E,  = 0,  (2.23) 

+ = (2.24) 


Boundary  conditions  for  scalar  Helmholtz  equations.  Let  H C denote 
the  transverse  cross-section  of  the  fiber,  and  50  the  core-cladding  interface  (50  is  a 
simple  closed  curve).  Since  0 is  translationally  invariant,  the  transmission  conditions 


(2.7-2.10)  can  be  expressed  as 

X n = E^^'^  X h on  dfl,  (2.25) 

7/(^1  X n = X n on  50,  (2.26) 

nj  E^^^  - n = nl  E^^'>  ■ h on  50,  (2.27) 

on  50.  (2.28) 


where  n — {nx,ny,0).  Then  from  (2.25-2.28)  the  z-components  E^  and  are  re- 
quired to  satisfy 

£;(1)  = on  50, 

i/i')  = on  50. 

These  conditions  are  not  enough  to  define  well-posed  problems  for  E^  and  i/^.  Addi- 
tional boundary  conditions  can  be  extracted  from  Maxwell’s  equations  and  equations 
(2.7-  2.10)  or  (2.25-2.28).  From  the  first  equation  in  (2.6),  we  find 

V X (1E(^^  - 1E(2))  X h = fu;/io(IH(^)  - IH^^^)  x h on  50 


0, 
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or 


(V  X X n = (V  X X n 


The  ^-component  of  this  equation  is 


JZ  tJJ-Jz 

-nx—^ n 


= —n 


This  implies  that 


dx  i/3{-nxEi^^  - UyEl^)) 

riy  - + ip{-nxEf'>  - UyE^^)  on  dO.. 


dn 


dn 


on  dVl. 


Similarly,  we  obtain 

[^V  X X M(2)]  X n = -iu;£o(E(i)  - E^^))  x h 

Ti\  Ho 

= 0 on  50. 


n? 


Hence,  ^ x h = ^ x h on  50.  The  2-component  of  this  equation  yields 

1 5i/i'>  1 


n 


j dn 


n\  dn 


on  50. 


So,  we  have  shown 


Proposition  2.3.3  The  weakly-guiding  approximation  of  Maxwell’s  equations  yields 
the  scalar  Helmholtz  or  2D  wave  equations  for  E^  and 


AE,  -f-  {k^n^  - ^^)E,  = 0 in 


(2.29) 


with  Ez  and  continuous  across  50  and 


AHz  + {Pn^  - f3^)Hz  = 0 in 


(2.30) 


with  Hz  and  continuous  across  50. 
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Remark:  The  boundary  conditions  on  and  are  slightly  different.  The  normal 
derivative  of  Ez  is  required  to  be  continuous  across  the  media  while  that  of  does 
not  have  be. 


CHAPTER  3 

PERTURBATION  OF  GUIDED  MODES  IN  CIRCULAR  FIBERS 


In  the  first  half  of  the  chapter,  approximations  of  eigenvalues  of  a family  of 
linear  self-adjoint  operators  using  the  classical  perturbation  theory  [20]  are  briefly 
considered.  These  approximation  methods  are  later  used  to  study  the  reduction  of 
the  signal  bandwidth  of  weakly-guiding  circular  fibers  whose  refractive  index  profiles 
are  slightly  deviated  from  the  optimal  truncated  parabolic  shape. 

3.1  Eigenvalue  Approximations 

In  this  section,  we  derive  some  useful  approximations  for  eigenvalues  of  a self- 
adjoint  operator  based  on  the  perturbation  theory.  Let  H he  a Hilbert  space  with  the 
inner  product  (u,v).  Let  L be  a self-adjoint  linear  operator  with  domain  of  definition 
D C H.  Let  G be  a domain  of  the  complex  e-plane. 

Definition  3.1.1  [20]  A family  L{t)  of  closed  operator  in  a Hilbert  space  H is  said 
to  be  analytic  of  type  (A)  over  G if  each  domain  of  definition  D{L{e))  = D is  inde- 
pendent of  e and  L{e)u  is  analytic  for  e € G for  every  u ^ D. 

Definition  3.1.2  A linear  operator  A is  said  to  be  relatively  bounded  with  respect  to 
L if  D(A)  = D[L)  = D and  there  exist  nonnegative  constants  a and  b such  that 

||Au||  <a||u||  + 6||Lu||  Vu  e D.  (3.1) 
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Let  Aq  6 IR  be  an  isolated  eigenvalue  of  L and  uq  be  the  corresponding  normalized 
eigenvector,  Le., 

Luq  = AoWo- 

If  |c|  is  sufficiently  small,  eigenvalues  A(e)  and  corresponding  normalized  eigenvectors 
u(c)  of  an  analytic  family  of  type  (yl),  {1(e)},  are  also  analytic  functions  in  e,  and 
moreover,  u(c)  and  A(e)  have  Taylor  expansions 

n(e)  = Uo  + + • ■ • , 

A(e)  = = Ao  + cA^^^  + + • • • , 

with  Z/(0)uo  = Luq.  Thus,  the  eigenvalue  equations 

L{t)u[t)  = A(e)u(c)  (3.2) 

can  be  written  as 

L{t){uo  + -)-•••)  = A(c)(uo  + + •••).  (3-3) 

Hence,  the  first-order  approximation  of  this  equation  is 

L[e)uo  = A(e)uo  -f  0(e).  (3-4) 

Thus,  we  obtain  the  first  simple  formula  which  is  widely  used  [28] 

A(e)  = Ao  -f  ([L(e)  — Zjuo,  uo).  (3-5) 

In  ([20],  p.421),  the  eigenvalue  problem  is  restricted  to  the  finite  system  of  eigenvalues 
- a finite  collection  of  eigenvalues  with  finite  multiplicities.  Hence,  the  perturbation 
theory  of  the  finite  dimensional  spaces  (for  matrices)  can  be  applied  without  much 
modification.  We  will  adapt  this  approach.  Let  A be  an  isolated  eigenvalue  of  the 
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unperturbed  self-adjoint  operator  L = L(0)  with  multiplicity  m.  Let  P be  the  asso- 
ciated eigenprojection  (here,  the  eigennilpotent  D is  identically  zero  because  of  the 
self-adjointness  of  L).  The  eigenvalue  A will  in  general  split  into  several  eigenvalues 
of  L{t)  for  small  e and  form  the  so-called  A-group.  Then  the  projection  P{(.)  for  this 
A-group  is  analytic  in  at  e = 0 G G.  Let  M(e)  :=  P{e)H  be  the  invariant  subspace 
under  the  perturbed  operator  L(t).  The  A-group  eigenvalues  of  L(e)  are  then  the 
eigenvalues  of  L(e)  in  the  subspace  M(e).  Thus,  in  order  to  determine  the  eigenval- 
ues in  the  A-group,  we  only  need  to  solve  the  eigenvalue  in  the  finite  dimensional 
subspace  M(e). 

Definition  3.1.3  The  weighted  mean  of  the  \-group  eigenvalues  of  L{e)  is  defined 
as 

A(e)  :=  Ur(L(,)P{c}) 
m 

= A + lir((i(e)-A)/>(£)). 

If  there  is  no  splitting  of  A so  that  the  X-group  consists  of  only  one  eigenvalue  A(e) 
with  the  same  multiplicity  m as  X,  we  have 

X{e)  = X{e). 

Proposition  3.1.4  Suppose  that  L has  a finite  number  of  isolated  eigenvalues 
{A,/ii,/i2, . . . ,/iiv},  where  the  eigenvalues  pk  are  different  from  X.  Let  P = P{0)  and 
{Pt}  be  the  corresponding  eigenprojections.  Then  the  weighted  average  eigenvalue 
A(e)  is  of  the  form 

A(e)  = A + eA(') -t- e^A^^)  + . . . 

where 

A(i)  = -trL(^)p{0), 
m 


(3.6) 


19 


and 


= -i-(r(iP>P(0)-Z,l‘>5i<‘IP(0)), 


Su=Y,-^ 


(3.7) 


fJ'k-  >< 


Proof:  Let  V be  the  total  eigenprojection  for  L,  i.e., 


LV= 

For  |c|  sufficiently  small,  we  have 

i(£)7’(e)  = Y. 

f‘i(Oe{A(e),Mlk(£)} 

where 

OO 

n=l 

Then  the  proposition  is  proved  using  the  standard  perturbation  theory  in  the  finite 
dimensional  case.  ■ 

Assume  that  the  eigenvalue  Aq  is  simple  (m  = 1).  Since  the  normalized  eigenvectors 
of  L form  an  orthonormal  basis  for  the  finite  dimensional  subspace  VH,  we  can 
express  the  formulas  in  the  previous  proposition  as 

A^i)  = (3.8) 

A»)  = (L(^'h.u)-  V ^ 

where  u and  vj  are  the  associated  orthonormal  eigenvectors  of  A and  /Xj,  respectively. 
In  fact,  the  operator  S is  of  the  form 


Sw  = ^ {fik  - Ao)  ^PjW 

/Aj5^Ao 


\/w  6 H. 
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The  size  of  |e|  determines  the  accuracy  of  the  approximations  in  (3. 5, 3. 8, 3. 9). 
We  will  use  these  formulas  to  compute  the  perturbed  eigenvalues  later.  We  now  derive 
an  upper  bound  for  |e|  such  that  the  formulas  in  (3. 5, 3. 8, 3. 9)  are  applicable.  The 
technique  used  here  is  based  on  the  theory  of  perturbation  of  a relatively  bounded 
operator. 

Definition  3.1.5  Let  L be  a closed  operator  in  H and  ( be  a complex  number.  If 
L — ( is  invertible  with 

R^L)  :=  {L  - C)-'  (3.10) 

then  ^ is  said  to  belong  to  the  resolvent  set  of  L.  The  operator-valued  fuction  R^{L) 
defined  on  the  resolvent  set  P{L)  is  called  the  resolent  of  L. 

Let  us  consider  a family  of  self-adjoint  perturbed  operators  L{t)  of  the  form 

L{e)  = L -I- cL^^\  D{L{e))  = D{L)  =:  D,  (3.11) 

where  is  a relatively  bounded  operator  with  respect  to  L, 

< a||u|| -t- 61|Lu||,  u ^ D,  (3.12) 

where  a and  b are  nonnegative  constants.  Let  ( be  a point  in  the  resolvent  set  P{L) 
of  L.  Since  L is  self-adjoint,  we  have  ([20],  p.272) 

ii^a^)ii  = \\{L-cr\\ 

= spr/?^(L) 

1 

Aecr(L)  K ~ 

1 


dist(C,cr(L)) 
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1 

~d' 

where  cr{L)  is  the  spectrum  of  L.  The  next  lemma  gives  a necessary  upper  bound  for 
|e|  so  that  a point  C in  the  resolvent  set  of  L is  also  an  element  in  the  resolvent  set 
of  the  perturbed  operator  T(c). 


Lemma  3.1.6  Let  L{t)  = L + be  self-adjoint  in  H and  satisfy  (3.12).  Let 
( be  a point  in  the  resolvent  set  of  L.  If 


a + 6(|C|  + d)  'sicr 


then  is  also  in  the  resolvent  set  of  L{e).  Moreover, 

OO 

RdL.)  ■.=  (i(0  - C)-‘  = Y, 

n=l 


where 


(3.13) 


(3.14) 


Bo  = 


(3.15) 


B„  = (-i)"flai)[i“'flf(i)r. 


(3.16) 


and 

l|fin||<^,  n = 1,2,3,...  (3.17) 

Proof:  Let  ( € P{L)-  the  resolvent  set  of  L.  In  order  for  ( € P(L(e)),  L(e)  — ( must 
be  invertible,  i.e.,  (L(e)  — exists  and  bounded.  We  write 

i(£)-C  = £i'»  + i-C 


= [/  + £iWfl£(i)](L-0. 


Since  6 P{L),  R<^{L)  = {L  — is  well-defined.  Hence  we  only  need  to  show  that 
I -t-  eL^^^R^(L)  is  invertible  for  some  c.  The  Neumann  series  of  [/  -f  cL^^^ R^{L)]~^ 
converges  absolutely  if 


|£|  ||i<‘>fl((i)||  < 1. 
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Since  is  L-bounded,  for  u E H,  we  have 


\\L^^'>Rc{L)u\\  < a\\R^{L)u\\+b\\LR^{L)u\ 


< (a||ii;c(L)||  + 6||I/i:^(L)||)  ||n| 


< 


U 


The  last  inequality  can  be  seen  from  noticing 


{L-OR^iL)  = I, 


so, 

LRc{L)  = CRdL)  + I,  (3.18) 


and,  hence 


||Li?c(L)||  < |C|Pc(i:)||  + 


a 


Therefore, 

||i'‘'flc(i)||  < i|a  + Kiel  + <!)]  =:  KO-  (3.19) 

Thus,  if  |e|  < the  operator  I + R(^{L)  is  invertible.  It  implies  that  L(e)  — C 

is  invertible  and  therefore  C is  also  a point  in  the  resolvent  set  of  L{e).  We  now  can 
express  (L{e)  — as  a Neumann  series: 

OO 

fl<(i(£))  = B<(i)  ^|il"H((L)]”£”  (3.20) 

n=0 

OO 

= fl((i)  + ^H<(i)l£'‘’fl<(i))"£”. 

n=l 


(3.21) 
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Let  us  denote 

Bq  — R({L), 

and 


Bn  = Rc{L)[L^^^Rc{L)]-,  n = l,2,.... 

Finally,  for  n = 1,2,..., 

p„ii  < iifl<(L)iiiii<‘ifl<(i)ir 

This  completes  the  proof.  ■ 

Definition  3.1.7  Let  Ao  G cr{L)  be  an  isolated  eigenvalue.  Let  A = <r(Z)  — Aq.  The 
isolation  distance  of  Aq  is  defined  as 


d — dist{Xo,  A). 


(3.22) 


Let  Ao  be  an  isolated  eigenvalue  of  multiplicity  m < oo  with  isolation  distance  d and 
uq  the  corresponding  normalized  eigenvector.  Let  F be  the  circle  center  Aq  and  radius 
|.  We  set 


(J:=m^^(C)  (3.23) 

where  ^(C)  is  defined  in  (3.19).  With  these  assumptions,  we  state  the  next  theo- 
rem which  can  be  used  to  estimate  the  coefficients  of  the  Taylor  expansions  of  the 
eigenvalues  A(e)  and  corresponding  normalized  eigenvectors  u(e). 


Theorem  3.1.8  Let  S be  as  in  (3.23).  If  |e|  < the  circle  F (center  Ao,  radius 
dj2)  encloses  exactly  m (repeated)  real  eigenvalues  A(e)  of  L{e)  and  no  other  points 
of  a{L{e)).  Moreover,  if  Ao  is  a simple  isolated  eigenvalue  (m  = 1),  then 


A(e)  = Ao  + eA*'*  + e^A*“*  + ■ ■ ■ , 


(3.24) 
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u(e)  = uo  + + • • • , (3.25) 

where 

n = l,2,...,  (3.26) 

||u(")||  < 2"(^",  n = l,2,....  (3.27) 

Proof:  Since  |e|  < ^ < any  point  on  the  circle  F is  in  the  resolvent  set  of  L{e)  by 
the  previous  lemma.  Then  from  (3.14), 

OO 

H<(t(e))  = ^B„£",  Bo  = BdL), 

n—0 

converges  uniformly  on  the  compact  set  F.  Hence,  the  eigenprojection  E{e)  can  be 
expressed  as 

E(e)  := 


1 

27r* 


Bnd( 


where 


Thus, 


— Eq  End^ , 

n—\ 


n 


1,2,.... 


\\E{e)-Eo\\ 


n=l 


1 °° 


n=l 


l._l , 1 c5  1 

2^1-St  ’ 21-eS^2' 


(3.28) 
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Since  E{c)  and  Eq  are  eigenprojections  with  ||£^(e)  — £^o||  < 1,  E{e){H)  and  Eq{H) 
are  isomorphic,  i.e.,  dim  E{e){H)  = dim  Eq{H).  Thus,  the  eigenspace  corresponding 
to  the  part  of  the  spectrum  of  L(c)  in  the  interior  of  P has  the  same  dimension  as 
the  eigenspace  of  Aq.  Hence,  if  Ao  is  a simple  isolated  eigenvalue(m  = 1),  there  exists 
exactly  one  A(c)  inside  P which  is  an  eigenvalue  of  L{e).  This  fact  is  important  in 
numerical  computation. 

Let 

E{e)uo 


u(e)  := 


ii^(o«oir 


Since  E{e)  is  an  eigenprojection. 


T(e)u(e)  = A(e)u(c). 


(3.29) 


Then, 


Eq  + En 


Uq 


u(e)  = 


n=l 


oo  \ \ 1/2 

£o+E  ^ I ^0  I 

Since  (EqUo^Uq)  = (uo,^^o)  = 1,  we  get 


u(e)  = 


EqUq  + ^ ^ e"£/„uo 

n=l 


1 + e"(£'„U0)  Uq) 


n=l 


1/2- 


(3.30) 


The  series  (3.30)  converges  for  |e|  < In  fact,  using  binomial  expansion  and 
we  see  that 


1 + ^ m- 


n=l 


E 

,m=0 


_ 1 
2 

m 


E I'l”'*” 


>vn=l 


1 — 


E 

m—O 


_ 1 
2 

m 


eS 


1 — eS 
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1 1 
“ ' (1  - rife)''' 

= (1  - • (1  - 2tS)-^'^  < (1  - 2e8)-\ 


which  is  striclty  less  than  1 provided  that  |c|  < ^-  Thus,  we  write 

OO 

u(e)  = uo  + ^ u„e" 

n=l 

where 


hn\\  < 2V. 


(3.31) 


This  shows  the  approximation  for  the  eigenvector  u(e).  Now,  we  derive  an  approxi- 
mation for  A(e).  Again,  since  L(e)u(e)  — X(e)u(e),  (L(e)  — Ao)u(e)  = (A(c)  — Ao)u(e). 
So, 

X ((L(e)-Ao)u(£),uo) 

A(C)  - Ao  — — 

(u(e),uo) 


((L(e)  - Ap).E(e)Mo,uo) 
(£;(6)uo,uo) 


(3.32) 


Since  L(e)  is  closed  and  Rc{L)  is  continuous  (bounded)  in  ( and  Ran(/2^(L(c)))  C 
D{L{e)),  we  have 

{L(€)-\o)E{k)  = (i(£)-A„)(-i/lJ((L(£))rfc) 


= LmdL{e))d<;  - A„  (-i  j(  Jk(LM)dcJ 

= ~£[/+wMem 


27 


-2^jf(C-Ao)/ic(i(£)MC 


17 

= --^  i(C-Ao)B„.iC 

S L 2>r./r  J 

oo 

=:  E'"C- 

n=0 

Since  Bq  — R(^{L),  Cq  = 0.  In  fact, 


2^  /[/  + iH((i)|<(C  - A„ 


— Z/JSo  ~ ^qEq 


= 0. 


For  n = 1, 2, 3, ... , 


l|C„||  < ^jf|C-Ac|||B„||MC| 
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I- 


(3.33) 


Finally,  from  (3.32),  we  obtain 

OO 

^ ^ £ (C*n^0j  ^o) 


A(e)  — Aq  — 


n=l 


1 + c”(£'„Uo,  Wo) 


n=l 


= ^e"(C„no,no)-^(-l)^ 


n=l 


m=0 


e'^{EnUo,uo) 


_n=l 


Since  IlCnll  < |(J",  II^Fnll  < right  hand  side  of  the  previous  equation  is  bounded 

by 


jOO  oo/oOc-nX”*  1 


cS 


n=l 


m=0  \n=l 


n=l 


2(l-c5) 


d 2{l -c5) 

- 2'2{l-lc5)  ^ 

j °° 

< 5E^"^” 

n=0 


Ee"i" 


(3.34) 


Hence,  from  (3.33)  and  (3.34),  we  have 

OO 

A(e)  = Ao  + y~^  A„e”, 


n=l 


where  |A„|  < This  completes  the  proof.  ■ 

In  Section  7,  these  approximations  are  applied  to  compute  the  perturbed  propagation 
constants  of  circular  fibers. 


3.2  Scalar  Wave  Equation 


In  Chapter  2,  we  derived  the  Helmholtz  equation  and  the  boundary  conditions 
(2.29)  for  the  longitudinal  component  EzB~'^^  of  the  time-harmonic  electric  field  E. 
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The  circular  symmetry  of  a fiber  allows  us  to  reduce  the  Helmholtz  equation  to 
an  ordinary  differential  equation  using  polar  coordinates.  Let  (r,9)  be  the  polar 
coordinates  in  Assume  that  the  index  profile  n is  radially  symmetric  and  2- 
invariant,  i.e.,  n = n(r).  Then  the  Helmholtz  equation  for  becomes 


dr 


de^ 


(3.35) 


Using  separation  of  variables,  we  write 


E,  = ^'(r) 


Then  ^'(r)  and  satisfy  the  equations 


W 


+ = 0, 


d^vlr  1 d<H 


m 


+ -^+  ^ = 0, 


(3.36) 

dr^  ' r dr  ' ^ 

where  m must  be  an  integer,  since  it  is  required  that  the  field  is  self-consistent  on 
each  rotation  of  $ through  27t.  Thus, 


_ / cos(a  -b  md), 

\ sin(a  -b  mO), 

where  a is  a constant.  The  differential  equation  (3.37)  for  the  function  ^(r)  is  called 
Sturm-Liouville  equation.  The  boundary  conditions  for  ^ are 


lim-\/r^’(r)  = 0, 

r->0 

lim  \/r^(r)  = 0. 

I — ►oo 

We  also  assume  that  ^'(r)  and  Vl’'(r)  are  continuous  in  (0,oo)  by  the  transmission 
conditions.  For  each  m,  the  propagation  constants  (3  of  guided  modes  have  to  satisfy 
the  double  inequalities  [32] 


krici  < (3  < kuco- 


(3.38) 
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This  gives  us  the  range  of  allowed  eigenvalues.  The  performance  of  a fiber  depends 
on  the  group  time  delays  for  the  propagating  modes. 


Definition  3.2.1  The  time  delay  per  unit  length  of  each  propagating  mode  is 


_ Id (3 
c dk 


(3.39) 


where  c is  the  speed  of  light  in  free  space,  and  k = 2n/X. 


The  time  delay  r can  also  be  defined  as  follows.  For  convenience,  (3.37)  is  rewritten 
as 

-i(rW')'  + - k^n^{r)  + ^ = 0.  (3.40) 

Let  (/?i,^'i)  and  (/?2,^2)  be  solutions  of  (3.40)  for  k = ki  and  k — k2,  respectively. 
Then,  by  integration  by  parts,  we  have 

/•OO  ^OO 

(/?2  — /3j ) / = (Aij  — fcj)  / n^'^i'^2f'dr. 

Jo  Jo 

So,  if  (3i  = /3{k),  /?2  = (3{k  + h),  and  h — > 0 and  since  the  eigenfunctions  ^(r,  fc) 
smoothly  depend  on  k,  (Chap.l,  Theorem  8.4  in  [7]),  we  have 

2/?-^  f '^^rdr  — 2k  [ n^(r)^^rdr. 
dk  Jo  Jo 


Hence,  the  time  delay  is 


r = 


\k  f^  n^(r)^^(r)rdr 


(3.41) 


c /3  '^>‘^(r)rdr 

With  this  integral  formula,  we  can  avoid  the  numerical  differentiation  in  calculating 


r. 


3.3  Perturbed  Fibers 


In  this  section,  a ’’perturbed”  index  profile  is  the  result  of  adding  a small 
perturbation  to  the  ’’unperturbed”  profile  whose  solutions  of  the  scalar  wave  equation 
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are  known.  If  the  deviation  from  the  known  profile  is  small  enough,  perturbation 

methods  provide  a good  approximation  to  the  solutions  of  the  perturbed  (unknown) 

fiber.  Define  the  perturbed  index  profile 

...  f n(r)  + eb(r)  0 < r < a, 
n(r)  = < 

{ rici  r > a 

where  b{r)  is  a bounded  function  and  c is  between  1%  and  10%  of  the  difference 
{rico  — Tici).  The  perturbed  eigenvalues  associated  to  the  perturbed  index  profiles 
can  be  approximated  by  the  formulas  (3.5),  (3.9),  and  (3.26)  where  L is  the  Sturm- 
Liouville  operator  defined  by 


m 

r2 


(3.42) 


and  the  perturbed  operator  is 


L{c)  L — e2n{r)b{r)  — e^b^{r) 


(3.43) 


Hence,  :=  2n(r)6(r)  and  :=  6^(r). 


3.4  Truncated  Parabolic  Profile 


The  truncated  parabolic  profile  is  given  as 


n^{R)  = 


nl^{l-2AR^)  0<R<1, 
nl,  R > 1, 


where  R = rfa  and  a is  the  core  radius.  The  scalar  wave  equation  is  normalized  to 

{ ^ }*(«)  = 0.  « e (»-  -)■ 

Let’s  define  the  following  parameters  [32]: 


V = 

(3.44) 

u = 

as/kV^-D\ 

(3.45) 

w = 

a^/?2-fc2n^. 

(3.46) 

32 


Substituting  these  parameters  in  the  above  equation,  we  get 


cP 


1 d 


m 


0 < /?  < 1, 


For  each  allowed  m,  the  (exact)  solutions,  normalized  at  i?  = 1,  are  [34] 


(3.47) 

(3.48) 


RM.AV) 

KmjWR) 

Km{W) 


0 < i?  < 1, 

R>\, 


(3.49) 


where  k = ^,  p = y,  and  / the  radial  index.  is  the  Whittaker  function  of 

the  first  kind  [Ij.  It  is  an  entire  function  in  z and  k.  Km{z)  is  the  modified  Bessel 
function  of  the  second  kind  [1].  To  find  the  propagation  constants  we  utilize  the 
continuity  conditions  on  '^{R)  and  '^'{R)  and  get  the  eigenvalue  equation 

(3.50) 


M'(V)  K’(W) 


Muv)  ^ " K^iwy 

Remarks:  The  Whittaker  function  of  the  first  kind  is  defined  as 
M«,^(z)  = iFt(^  + /u  - K,2fx  + l,z) 


= e 


E 


n=0  ”• 

where  (p)„  = p(p  + l)(p  + 2) . . . (p  + n — 1).  This  power  series  converges  for  all 
finite  values  of  z and  k [33].  Hence,  this  function  can  be  easily  computed  from  this 
definition. 

If  the  time  delay  r is  computed  using  the  integral  formula  (3.41),  it  is  convenient  to 
use  the  integral  identities  of  K^(W R)  [32] 


/: 


i 


OO  1 

Kl(WR)IUR  = -(/f?(W')  - K^iW)), 
Kl{WR)RiR  = - Kl(W)). 
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3.5  Distortions 

In  this  section,  we  consider  several  different  types  of  distortions  of  the  trun- 
cated parabolic  profile.  All  of  the  following  distortions  are  simulated  by  adding  a 
perturbation  function  p{R)  to  the  truncated-parabolic  profile  n{R). 

(i)  The  ripple  distortions  (without  the  exponential  term)  is  simulated  by 

p{R)  =tsin{2nNR).  (3.51) 

The  parameters  are  the  amplitude  e which  is  between  1%  and  15%  of  the  difference 
«co  — ^ci  and  the  number  of  ripples  which  ranges  from  1 to  20. 

(ii)  The  ripple  distortions  with  the  exponential  term  decay  exponentially  away  from 
the  fiber  axis  and  are  simulated  by 

p{R)  = esm{2nNR)exp{—4R^).  (3.52) 


(Hi)  The  dip  distortion  often  occurs  when  fibers  are  made  by  the  MCVD  process. 
This  disortiion  can  be  simulated  by 


P{R)  = 


nco(l  — 2A/?^)2  — 6 • exp 


Rl 


(/?o  - Ry 


+ 1 0<R<Ro, 


(3.53) 


0 elsewhere 

where  b is  between  5%  and  10%  of  rico  — rid  and  Rq  is  between  5%  and  10%  of  the 
core  radius. 

(iv)  The  bulge  distortions: 

0.001  sin{57r(/2  — /?o)}  Ro  < R < Ro  + 0.2, 

0 elsewhere 


p(R)  = 


(3.54) 


where  the  parameter  Rq  is  the  location  of  the  bulge. 
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3.6  Frequency  Response  and  Bandwidth 

There  are  at  most  finite  number  of  propagating  modes  in  a fiber.  Each  mode 
travels  with  different  speed,  hence  arrives  at  the  end  of  a Ifcm-long  fiber  at  a different 
time.  In  [24]  and  [25],  the  impulse  response  is  approximated  as  follows.  The  interval 
of  the  arrival  times  of  all  propagating  modes  (time  spread)  is  divided  into  n time  slots, 
for  some  suitable  n,  the  number  of  guided  modes  arriving  in  each  time  slot  is  counted. 
The  result  is  a step-function  approximating  the  impulse  response.  The  frequency 
response  of  a fiber  is  defined  as  the  absolute  modulus  of  the  Fourier  transform  of  the 
impulse  response  step-function.  Then,  the  signal  bandwidth  is  the  half-maximum 
frequency  of  the  frequency  response.  However,  if  the  time  spread  is  very  large,  this 
approximation  of  the  impulse  response  yields  inaccurate  bandwidth.  An  alternative 
approach  is  given  in  [9].  Each  mode  is  treated  as  an  independent  delta  function.  Thus, 
the  pulse  at  1km  is  the  distribution  of  delta  functions  (spikes)  with  appropiately 
weighted  amplitudes.  More  precisely,  if  r,,  i = 1, 2, . . . , A^,  is  the  arrival  time  of  each 
propagating  mode  and  Wi  is  the  corresponding  modal  weight,  the  pulse  T is  described 

by 

N 

m = E Wi5(t  - n)  (3.55) 

i=l 

where 

The  frequency  response  F{ui)  is  the  absolute  value  of  the  Fourier  transform  of  the 
pulse  r. 


F{co) 


/oo 

T{t) 

•OO 

N 


dt 


_1_ 

2^ 


,g27riTi<ii 


1=1 


(3.56) 
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The  fiber  bandwidth  is  defined  the  same. 

The  truncated-parabolic  fiber  used  here  is  described  by 

a = 31.25/zm, 

A = O.Sbfim, 
rico  = 1.48/im, 
rict  = lA6fj,m. 

This  particular  fiber  supports  210  guided  modes,  whose  signal  bandwidth  is  1413MHz- 
km.  All  of  the  slightly  perturbed  fibers  considered  here  also  support  the  same  num- 
ber of  guided  modes  but  their  signal  bandwidths  are  reduced  very  significantly.  This 
can  be  seen  in  the  numerical  experiments  in  the  next  section. 

3.7  Numerical  Experiments 
3.7.1  First  Order  and  Second  Order  Approximations 

We  first  compare  the  eigenvalues  found  by  the  first  and  second  order  using 
(3.5)  and  (3.9)  with  those  found  by  the  finite  element  method  (discussed  in  the 
next  chapter).  Here  we  use  the  ripple  distortion  without  the  exponential  term  with 
e = 1%  • {rico  — riel)  = 2 • 10“^  and  1 ripple  and  the  operator  L{e)  = L -f  -(- 
defined  in  (3.43).  For  the  finite  element  error  in  eigenvalues  to  be  less  than  10~^°, 
it  requires  1024  mesh  points  for  the  fundamental  mode  (m  = 0,/  = 0)  and  8192 
mesh  points  for  the  higher  mode  (m  = 0,/  = 6),  while  the  formulas  (3.5)  and  (3.9) 
are  easily  computed.  Therefore,  the  perturbation  methods  are  much  faster  and  less 
expensive  especially  for  smaller  propagation  constants  (*.e.,  closer  to  the  cut-off). 

For  the  application  of  (3.26),  we  rewrite  the  operator  as  = ‘ln{r)h{r)-\- 
eb^{r).  For  m = 0,  the  isolation  distance  of  the  largest  exact  eigenvalue  of  L,  Aq  = 
10.9348901,  is  d = 0.03.  We  have  “f*  ^||^||  where  a = 2rico  + e,  since 
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Table  3.1:  Eigenvalues  by  First  and  Second  Order  and  FEM  for  Ripple  Distortion 


Unperturbed 

First  order 

Second  order 

FEM 

10.9348901 

10.9358125 

10.9358229 

10.9358230 

10.8717622 

10.8717934 

10.8718034 

10.8718035 

|^(^)|  ^ 1 ^•nd  6 = 0.  So, 


(5  = 


a 

d 


2- 1.48  + .0002 
d 


296.02. 


Therefore,  the  perturbed  eigenvalue  A(e)  corresponding  to  Aq  can  be  approximated 
by  the  Taylor  series 

A(6)  = Ao  + eA(i)  + c2a(2)  + . . . 

I 

with  c = 2 • 10-''  and  |A(")|  < = 5 • 10“^  x 296".  Hence, 


|A(e)-Ao|  < 9 • 10-'' + 0(e2). 


This  bound  agrees  with  the  result  in  the  above  table. 
3.7.2  Bandwidths  of  Perturbed  Fibers 


Using  the  perturbation  method,  we  compute  the  perturbed  eigenvalues.  Thus, 
the  bandwidths  of  perturbed  fibers  can  be  approximated.  We  assume  that  all  modes 
carry  the  same  modal  weight 


Wi  = 


1 

~N 


where  N is  the  total  number  of  propagating  modes.  As  we  can  see  from  the  tables, 
fibers  suffer  a serious  reduction  in  bandwidths  when  their  index  profiles  are  slightly 
deviated  from  the  optimal  shape.  For  example,  in  the  table  for  ripple  distortions 
without  the  exponential  term,  the  bandwidth  of  the  truncated  parabolic  fiber  is 
dropped  from  1413  MHz  • km  to  449  MHz  ■ km  when  its  profile  is  distorted  with  only 


37 


one  ripple  of  amplitude  e = 1%  • {rico  — nd)  = .0002.  We  plot  the  index  profiles  of 
different  perturbed  fibers  and  the  associated  impulse  responses. 

Table  3.2:  Bandwidths  of  Ripple  Distortions  without  exponential  term:  p(R)  = 
esm{2TrN  R). 


t 

N BD{MHz  ■ km) 

1%  • («co  - «d)  1 449 

1%  • 

5 

174 

1%  * {tIqq  ~ 

10 

194 

1%  • {rico  - rici) 

15 

1119 

1%  • {rico  - rid)  20 

1437  (noise) 

Table  3.3:  Bandwidths  of  Ripple  Distortions  with  exponential  term:  p{R)  — 

esin{2n NR)  exp(  — -^). 


e N BD{MHz  ■ km) 

1%  • {rico  - rid)  10 

557 

2%  • {rico  - rid)  10 

298 

5%  • {rico  - rid)  10 

121 

10%  • {rico  - rid)  10 

62 
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T runcat^d  parabolic  profila 


Figure  3.1:  Truncated  parabolic  profile. 


Figure  3.2:  Frequency  reponse  of  the  truncated  parabolic  profile. 
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Ripple  ,opsilon*1  %,N-5 


Figure  3.3:  Ripple  distortion  profile,  = 5,  e = 1%  • {rico  — nd). 


Figure  3.4:  Frequency  responses  of  the  optimal  profile  (solid)  and  of  the  ripple  dis- 
tortion (dotted)  with  N = 5,  e = 1%  • {uco  — rici)- 
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Ripple  with  exponential  term,epsilon-5%,N«10 


Figure  3.5:  Ripple  distortion  with  exponential  decay,  N = 10,  e = 5%  • {uco  — nd). 


Figure  3.6:  Frequency  responses  of  the  optimal  profile  (solid)  and  of  the  ripple  dis 
tortion  (dotted)  with  exponential  decay,  = 10,  e = 5%  • {uco  — rici). 
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Dip.  R0-.05,  b-5% 


Figure  3.7:  Dip  distortion  profile,  Rq  — .05,  e = 5%{rico  — nd). 


Figure  3.8:  Frequency  responses  of  the  optimal  profile  (solid)  and  of  the  dip  distortion, 
Ro  = .05,  b = 5%{tIco  - rici). 
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Bulge,  RO-0 


Figure  3.9:  Bulge  distortion,  i?o  = 0- 


Figure  3.10:  Frequency  response  of  the  optimal  profile  (solid)  and  of  the  bulge  dis- 
tortion (dotted)  , = 0. 
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Table  3.4:  Bandwidths  of  Dip  Distortions:  p{R)  = nco(l  — 2A/2^)2  — b • 

exp  + l)  , 0 < i?  < i?o. 


Table  3.5:  Bandwidths  of  Bulge  Distortions:  p{R)  = esin{57r(/?—  Rq)},  Ro  R 
Rq  + 0.2. 


CHAPTER  4 

A FINITE  ELEMENT  METHOD  FOR  CIRCULAR  FIBERS 


The  propagating  fields  in  a weakly  guiding  circular  fiber  are  governed  by  the 
scalar  wave  equation  which  is  a singular  Sturm-Liouville  eigenvalue  problem  on  the 
unbounded  interval  (0,oo).  The  finite  element  method  developed  by  Hohn  [17],  is 
generalized  to  solve  the  scalar  wave  equations. 

4.1  Unbounded  Domain  or  Non-Compact  Problem 


We  recall  that  Maxwell’s  wave  equations  of  electromagnetic  fields  in  a weakly 
guiding  circular  fiber  can  be  approximated  by  scalar  wave  equations  by  ignoring  the 
polarization  effects.  We  again  assume  that  the  index  profile  n is  radially  symmetric 
and  longitudinal  invariant,  i.e., 


n 


n^(r) 


n: 


0 < r < a, 
r > a 


(4.1) 


where  a is  the  core  radius,  n{r)  is  real-valued,  bounded,  piecewise  continuous,  and  rid 

(a  positive  constant)  is  the  refractive  index  in  the  cladding.  Therefore,  the  cylindrical 

coordinates  (r,  6,  z)  are  used  where  z is  the  fiber  axis.  The  longitudinal  component 

of  TE  fields,  denoted  here  as  Em,  can  be  expressed  as 

p - / ^m(r)cos(m6>)exp(-z7?z)  even  mode 

\ ’^m{r)s\n{m6)exp{—i/3z)  odd  mode  ^ ' 

where  /?  is  the  propagation  constant,  and  is  the  solution  of  the  scalar  wave 

equation 


u"  -I-  — -I-  ^ k^n^{r)  — ^ u = r G (0,  oo) 


m 


(4.3) 
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where  k = 27t/A  is  the  wave  number.  This  second-order  differential  equation  is 
classically  known  as  the  singular  Sturm-Liouville  eigenvalue  problem  because  of  the 
unbounded  domain  and  singularities  at  r = 0 and  r = oo.  In  the  traditional  fi- 
nite element  methods,  the  unbounded  domain  is  truncated  to  a finite  domain,  i.e., 
u{r)  = 0,  Vr  > To  > 0,  and  the  finite  dimensional  spaces  14  are  used  to  approximate 
the  eigenfunctions.  This  truncation  method  works  well  provided  that  the  eigenfunc- 
tions decay  rapidly  like  e""'"  for  a > 0 sufficiently  large.  For  a > 0 small,  the  finite 
dimensional  approximate  spaces  might  perform  poorly.  The  second  disadvantage  of 
this  method  is  that  if  the  original  problem  is  of  the  non-compact  case  (e.g.  Rellich- 
Garding  compact  embedding  condition  [3],  [19]  does  not  hold),  the  approximating 
eigenvalues  Xh  of  the  discrete  problem  (compact)  might  converge  to  a non-eigenvalue 
of  the  continuous  problem  [31].  This  phenomenon  is  known  as  the  existence  of  spu- 
rious solutions.  Such  non-compact  differential  operators  on  a bounded  interval  has 
been  studied  in  [27]  and  on  an  unbounded  interval  in  [17].  We  now  briefly  describe 
the  methods  developed  in  [27]  and  [17].  Let  77  be  a Hilbert  space  and  V C H a 
closed  subspace.  Consider  the  generalized  eigenvalue  equation 

Lu  — XMu  (4-4) 

where  [L,D{L))  and  [M,D{M))  are  differential  operators  with  D[L)  C V and 
D{M)  C H.  Let  (Lh,D{Lh))  be  the  corresponding  differential  operator  for  the 
discrete  problem.  Assume  L and  Lk  are  invertible  with  the  inverses  L~^  and 
respectively.  Denote  T L~^M  and  Th  :=  Mh-  So  the  general  eigenvalue  prob- 
lem is  reduced  to  find  eigenvalues  of  the  bounded  operators  T and  Th-  If  the  finite 
element  subspaces  are  finite  dimensional,  Th  is  finite  rank,  hence  compact.  The  usual 
a.ssumption  of  the  norm  (uniform)  convergence  [12]  jjr  — 24][  — 0 as  — >•  0 excludes 
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the  situation  where  T is  non-compact.  In  [27], the  eigenvalue  problem  is  defined  on 
a bounded  interval  and  a weaker  condition  of  strong  convergence  \\Tu  — Thu\\  0 
instead  of  the  norm  convergence  is  assumed.  As  in  the  compact  case,  each  isolated 
eigenvalue  of  T can  be  approximated  by  a sequence  of  eigenvalues  of  the  compact 
operators  Th,;  however, it  is  not  true  in  general  that  every  sequence  of  approximate 
eigenvalues  will  converge  to  an  (isolated)  eigenvalue  of  T.  In  [17],  the  eigenvalue  prob- 
lem is  defined  on  an  unbounded  interval  with  singularities  at  the  end-points.  The 
approximate  subspaces  14  are  chosen  to  be  infinite  dimensional  as  follows.  Let  / be 
a bounded  interval.  On  7,  14  consists  of  the  usual  finite  element  functions.  Outside 
of  7,  14  is  identical  to  V.  This  method  eliminates  the  two  disadvantages  mentioned 
above.  However,  it  requires  the  fundamental  set  of  solution(s)  at  each  end-point  to 
be  known.  In  this  chapter,  we  extend  the  result  of  this  work  to  solve  our  scalar  wave 
equation  on  (0,  oo)  where  only  the  fundamental  solution  at  infinity  is  available. 

4.2  Variational  Formulation  for  E,  Field 
Let  77  :=  7/^(0,  oo,r)  be  the  Hilbert  space  equipped  with  the  inner  product 

/•CO 

{u^v)  :=  {u,v)h '■=  I uvrdr,  (4-5) 

Jo 

and  hence  77-norm  is  defined  as 


Consider  the  differential  operator  L 

7/U  :=  — -(ru')' -t- — A:^n^(r)|  u (4.7) 

with  the  domain  D(L)  = {u  G CHO,  oo)  : Lu,  u,  u'  G 77,  and  ^ ^ H if  m ^ 0}.  Then 
the  scalar  wave  equation  can  be  expressed  in  terms  of  L as  {L  + (J^)u  = 0. 
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Consider  the  following  variational  formulation 


u'v'rdr  + j ^ — k^n^{r)  > uvrdr  = —(3'^  J uvrdr.  (4.8) 


Let 


V :=  {u  E H : u'  E H,  lim\/ru(r)  = lim  \/ru(r)  = 0}. 

r->0  ‘ 


r—^oo 


(4.9) 


V is  a closed  subspace  of  H,  V = H,  and  is  equipped  with  the  weighted  norm 

||u||v.=  (r(|«P+|uf)rdr)’.  (4.10) 

Proposition  4.2.1  If  u ^ V , — £ H . 

r 

Proof:  First,  we  show  the  proposition  for  u G C^(0,  oo)  D V,  then  by  the  density  of 
in  V we  are  done.  Since  limr_4.oo  u{r)  = 0,  we  have 

u{r)  1 


1 

= --  J u'dt. 


Let 


1 7°° 

g{r)  - I u'{t)y/tdt. 
r Jr 

Then  by  Hardy-Littlewood’s  inequality  [10],  we  have 

/•CO  poo 

I 9^dr  < 2 / u'^rdr. 
Jo  Jo 

Using  this  inequality,  we  get 


L = / H/ 

2 


- / M/ 

<-a: 


u^rdr. 


(4.11) 
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'LL 

Hence,  — < 2||u'||  < oo.  ■ 
r 


Remark:  The  inequality  (4.11)  complements  the  well-known  Hardy’s  inequality 


Therefore,  a{u,  v)  is  a bilinear  form  defined  on  V x V for  m = 0, 1, 2, . . . 

Proposition  4.2.2  The  bilinear  form  a{u,v)  is  symmetric,  continuous  on  V x V 
and  bounded  from  below  by  — 

Proof:  Let  u,v  E.V.  It  is  clear  that  a(-,-)  is  symmetric.  By  Holder’s  inequality  we 
have 


So,  the  bilinear  form  a(-,  •)  is  symmetric,  and  bounded  from  below  by  the  constant 


Remark:  If  q > the  bilinear  form 


< ||u'|||H|-Hm^Hr||||u/r|KA:^||n^|U||uliy^ 


< C'(m)||u||v||n||v. 


Hence,  a{u,v)  is  continuous.  Also, 


b{u,v)  :=  a{u,v)  -f  a{u,v) 


(4.13) 


is  continuous  and  V-elliptic.  The  space  V with  the  inner  product  b{u,  v)  is  a Hilbert 


space. 
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By  a representation  theorem  ([20],  theorem  VI. 2. 2. 6),  there  exists  a unique  self-adjoint 
linear  operator  A with  B(A)  C D{a)  = V and 

a(u, u)  = (Au, u)  yuED{A),v^V  and  (4-14) 

(^u,u)  > (4.15) 

Proposition  4.2.3  The  spectrum  of  A,  cr{A),  is  a subset  of  the  interval  [— oo) 
The  essential  spectrum  of  A,  aess{A)  is  the  subinterval  [—k^nh,oo). 

Proof:  Since  A > a{A)  C [-Pn^„^,oo). 

By  definition,  A = L on  D{L).  By  the  transformation  w{r)  = y/ru{r),  we  obtain  the 
differential  equation 

/ — — \ 

—w"  -f-  I ^ — k^n^{r)  j w = (4.16) 

We  have  limr^oo  ~ k'^n^fr)  — —k’^nh.  Hence,  (Tess{A)  — [— A;^n^;,oo)  [29].  ■ 

Remark:  The  propagation  constants  /3  of  propagating  modes  must  satisfy 

kn^i  ^ ^ kn^Yifix' 

For  Q > the  non- homogeneous  equation 

[A  + a)u  = f (4.17) 

has  a unique  solution  u € D{A)  for  any  f ^ H,  because  —a  is  in  the  resolvent  set 
p{A).  Hence,  the  resolvent  of  yl  -(-  a exists  and  is  bounded,  denoted  by 

T:=  R^{A)  :=  (T  + a)"\ 

The  operator  T is  defined  everywhere  on  H with  the  range  D{A).  If  p is  an  eigenvalue 
of  A + a,  A :=  4 — Q,  is  an  eigenvalue  of  T.  Hence,  the  problem  of  finding  eigenvalues  of 
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the  differential  operator  A + a reduces  to  finding  eigenvalues  of  the  bounded  operator 
(Green  operator)  T.  The  corresponding  eigenfunctions  of  A + a and  of  T are  identical. 
Embeddings  in  weighted  Sobolev  spaces.  The  Hilbert  space  H = L^(0,oo,r)  is 
an  example  of  a power-type  weighted  Sobolev  space  on  (0,  oo)  [21].  Define 

oo,  r)  = € L^(0,  oo),  j = 0, 1, . . . , m}.  (4-18) 

Lemma  4.2.4  [SlJ  The  set  A = {4>  ^ C'°°(0,oo)  : lim  u{r)  = 0}  is  dense  in 

I — ^oo 

L^(0,  oo,  r). 

Corollary  4.2.5  The  set  B = {4>  ^ C°°(0,  oo)  : lim  y/ru(r)  = 0}  is  dense  in 

r— foo 

T^(0,  oo,  r). 


Lemma  4.2.6  [2]  If  u{r)  € C^(0,d),  d < oo  and  u'  G L^(0,d,r),  then 

lim  \/r  u(r)  = 0. 

r->^0 

Corollary  4.2.7  If  u{r)  G C^(0,oo)  and  u'  G L^(0,oo,r),  then 

lim  yJr  u(r)  = 0. 

r^O  ' ' 


Lemma  4.2.8  Let  u{r)  G B where  B is  defined  in  corollary  4-2.5.  Then 


poo  pd 

I \u\^dr  ^ 2 / |u||u'|rdr. 

Jo  Jo 


Proof:  By  integration  by  parts  and  corollary  4.2.7,  we  have 


poo  pc 

I \u\^dr  = —2 
Jo  Jo 


lu|  — |u(r)|rdr. 
dr 


Thus, 


poo  poo 

I < 2 / |u||u'|rJr. 

Jo  Jo 


(4.19) 


This  completes  the  proof.  ■ 
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Corollary  4.2.9  If  u E //^(0,oo,r)  then  u € L^{0,oo)  and 

lhlU2(0,oo)  < 2||u'||£,2(0,oo,r)-  (4.20) 

Proof:  The  proof  follows  directly  from  the  previous  lemma.  ■ 

Finally,  for  completeness,  we  state  a lemma  in  [2]. 

Lemma  4.2.10  Let  u G C^(0,d).  If  p>\,  u > 0,  then 

o pd  pd 

sup  |u(r)|^  < - I \u\^dr I 

(o,d)  d Jq  Jq 

sup r"|u(r)|P  < — + 2p  f |u|^~^|u'|r''dr. 

(0,d)  d Jq  Jq 

Corollary  4.2.11  If  u G C^(0,oo),  then 

poo 

sup  |u(r)|^  < / |u||u'|dr, 

(0,oo)  Jo 

roo 

sup  |u(r)|^r  < 4 / luHu'Irdr. 

(0,oo)  J 0 

4.3  Finite  Element  Approximations 

Let  14  be  a subspace  of  V,  not  necessarily  finite  dimensional.  Let’s  denote  Hh 
as  the  completion  of  14  in  H.  By  the  representation  theorem,  for  a{uh,  Vh)  defined 
on  14  X 14,  there  exists  a unique  self-adjoint  operator  Ah  with  D{Ah)  C 14  which  is 
bounded  from  below  by  and  satisfies 

(AhUh,Vh)  = a{uh,vh),  yuh  e D{Ah),Vh  e Vh.  (4.21) 

Ah  is  not  necessarily  compact  since  14  is  not  finite  dimensional.  The  inverse  of  Ah  + a 
exists  and  bounded  in  Hh  for  a > denoted  by  24.  This  implies  that  for  every 

/ G Hh,  the  unique  solution  of  the  non-homogeneous  equation 


{Ah  + a)u  - f 


(4.22) 
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is  Uh  :=  Thf.  Equivalently,  the  function  satisfies  the  variational  equation 

a{uh,Vh)  + a{uh,Vh)  = if.Vh),  ^Vh  e Vh-  (4.23) 

Remarks:  There  are  three  projections  of  V onto  Vh  that  are  useful  for  our  purposes. 

1.  The  orthogonal  projection  Ph  : H Vh  defined  as  follows.  For  each  f E Hh, 

{Phf,Vh)  = {f,Vh),  Vu/,  e 14.  (4.24) 

With  this  projection,  we  can  extend  the  operator  Th  to  H by 

Thf  :=  ThPhT,  V/  e H.  (4.25) 

In  fact,  the  solution  Uh  :=  Thf  of  the  variational  equation  (4.23)  remains  the  same 
for  every  f ^ H. 

2.  The  interpolating  projection  Wh  : H Sh{I)  where  Sh(I)  is  a finite  dimensional 
(closed)  subspace  of  Vh  consisting  of  polynomials  defined  on  a bounded  interval  I. 
We  will  come  back  to  Sh  later. 

3.  The  eigenspace  projection  Q : H A^/(Ao)  C V. 

Let  Ni{Xo)  be  the  subspace  of  generalized  eigenfunctions  of  T associated  with  the 
eigenvalue  Aq  of  (algebraic)  multiplicity  /.  Let  F be  a simple  closed  curve  about  Aq 
such  that  the  region  bounded  by  F and  containing  Aq  intersects  with  the  spectrum 
of  T only  at  Aq.  The  projection  Q is  defined  by  the  Cauchy-type  integral  [20] 

Q ■■=  f (T  - (4.26) 

Hence,  if  Aq  is  an  eigenvalue  of  T,  TQu  = XqQu,  Vu  € H.  Similarly,  we  can  also 
define  the  projection  Qh  • H Ni{Xq)  where  Aq  is  an  eigenvalue  of  the  operator  Th- 
We  hope  that  QhU  — > Qu  as  /i  — > 0. 
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The  subspaces  Vh  are  assumed  to  satisfy  the  following  approximation  properties. 
Approximation  Properties  of  Vh 

A1.||T  — T/i||  < where  lim/i_4.o  7/i  = 0 (i.e.,  uniform  convergence). 

A2.  Let  Ao  7^  0 be  an  isolated  eigenvalue  of  T of  multiplicity  /.  Define  the  gap 
between  the  eigenspace  associated  with  Aq  and  Vh  as  [20] 


e{h)  :=  sup  dist(u,V/,) 

u€JV((Ao) 

IMI<i 

:=  sup  inf  ||u  — u;,||.  (4.27) 

“€N,(Ao) 

IMI<i 

In  particular,  if  Aq  is  a simple  eigenvalue  of  T and  uq  is  the  corresponding  eigenfunc- 
tion, then  e{h)  = dist(uo)  Vh). 

We  assume  that 


lim  e(h)  = 0.  (4.28) 

h-*0 

The  assumption  (Al)  means  that  for  any  f E H the  solutions  u :=  T f and  Uh  Thf 
of  the  non-homogeneous  problems  ((A-(-a)u,u)  = (/,  u)  and  {(Ah  + oi)uh,Vh)  = {f^Vh) 
respectively,  satisfy 

< 7/i-  (4.29) 


In  the  next  lemma,  we  consider  the  regularity  of  solutions  u 6 D{A)  of  the  equation 
Au  = /,  i.e., 


—u"  -- — ( ^ \ u = f,  0 < r < oo. 


Lemma  4.3.1  Let  f £ H . If  u G D[A)  satisfies  Au  = f,  then  u G /7^(0,oo,r). 


Proof:  Let  u G D{A)  be  a solution  of 


-u I k n — 


r 


u = f,  0 < r < oo, 
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or,  equivalently, 


m 


{ru')  — ( k^rv — 1 u • r = f • r. 


Since  limr  • u'{r)  — 0,  we  have 


r->0 


-u'  = -J  ^k^n^(s)  — u{s)s  ds  + - J f{s)sds. 


Since  lim\/r  • u(r)  = 0,  we  further  have 


r-fO 


U 


- J ( A;^n^(s) —J  u{s)sds  (/<  + y - J f{s)sds  dt. 


u 


By  differetiating  the  integral  representation  of  u twice,  we  obtain 

" = ^ J (k'^n'^is)  - u{s)sds  - (^k^n'^(s)  - u{r)  + ^ J f{s)sds  - f{r). 

Using  Hardy-Littlewood’s  inequality,  z.e., 

1 rx 


f{s)ds 


L2(0,oo) 


< 2||/||l2(o.oo), 


we  obtain 


< ||-  / ( k^n^{s)  — ^ ) u{s)\l-y/sds 


+ 11“/  fi^)\l-Vsds 


L2 


+ 


L2 


m 


k'^n^ ) u^/r 


L2 


+ ll/^ll 


L2 


— 11“  / ( u{s)-^di 

Ik  ^0  V ^ J 


+ ||-  [ f(s)y/^ds 


Z.2  II  ' ^0 


+ 


L2 


7 2 2 ^ \ r 

k n — ) uy/r 


L2 


+ ll/^ll 


L2 


< Cdl^V^IlL-  + Il4v^||z,>  + ||/V?IU.)- 


Since  u G D{A),  ||>/ru"|k2  < oo.  ■ 
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4.4  Convergence  of  Eigenvalues  and  Eigenfunctions 

First,  we  restate  the  convergence  results  of  the  approximate  eigenvalues  and 
their  corresponding  eigenfunctions.  From  now  on  we  assume  that  ||T  — 7^11  — > 0 as 
h ^ 0. 

Lemma  4.4.1  Let  Aq  be  an  isolated  eigenvalue  of  the  self-adjoint  operator  T with 
multiplicity  1.  Then  for  each  e > 0 there  exists  ho{e)  > 0 such  that  for  any  h < ho(c) 
the  spectrum  of  the  self-adjoint  Th  consists  of  I eigenvalues  A)j  (counting  multiplicity) 
and  they  satisfy 

|A*,-A|<c,  i = l,2,... 

Proof:  This  is  a well-known  result  in  ([20],  1V§3.4).  ■ 

Lemma  4.4.2  Let  be  a sequence  of  isolated  eigenvalues  of  Th  converging  to  A. 
Then  A is  either  an  isolated  eigenvalue  ofT  or  —k^nh. 

Proof:  Since  T and  Th  are  self-adjoint  bounded  operators,  by  ([20],  V §4,  theorem 
4.10) 

sup  dist(A/i,c7(r))  < ]|r  - T/i]], 

A/,6<r(Th) 

sup  dist{X,  a{Th))  < ||r  - Th\\. 

Aecr(T) 

So,  if  A were  in  the  resolvent  set  of  T,  there  exist  7 > 0 such  that  dist(A,  <r(r))  > 
7.  However,  jjT  — Th\\  0,  there  is  ho  > 0 such  that  dist(A;jo, (r(r))  < 7/4  and 
sup;j<;,g  \Xh  — A|  < 7/4.  Thus,  it  is  a contradiction.  If  A is  a point  in  the  essential 
spectrum  of  T,  it  has  to  be  the  point  —k^n^  because  cress(^)  = °°)-  • 

Lemma  4.4.3  Let  A be  an  isolated  eigenvalue  of  T of  multiplicity  1.  Let  {A/,}  be 
a sequence  of  eigenvalues  of  Th  such  that  lim/i_>.o  Xh  = X and  {uh}  C Vh  be  the 


56 


sequence  of  the  associated  eigenfunctions.  Then  there  exists  a subsequence  of  {u/,} 
that  converges  to  u E V,  u ^ 0,  and  u is  an  eigenfunction  of  T . 

Proof:  Let  P be  a simple  closed  curve  containing  A in  its  interior  and  no  other  points 
of  the  spectrum  of  T.  Let  D be  the  region  enclosed  by  P.  For  h small  enough,  there 
are  / eigenvalues  of  T/,  (counting  multiplicity)  such  that  they  are  all  in  D.  Hence, 
the  difference  of  the  eigenprojections  Q and  Qh  can  be  written  as 

Q-Q,  = -^J^(RdT)-Rc(TkM.  (4.30) 

Therefore, 

IIQ-O-II  < f;  £\\R({T)  - 

Since  the  curve  P is  compact  and  R^{T),  R^{Th)  are  continuous  functions  on  the 
resolvent  sets  of  T and  Th,  respectively,  we  conclude  that  HQ— Qa||  0 for  ||r— T;i|| 

0.  Hence,  ||Qu  — Q/iu||  — > 0,  Vu  E H.  Therefore,  there  is  a subsequence  in  {Qhu}, 
denoted  by  Uh,  that  converges  to  an  eigenfunction  u E {Qu}.  ■ 

Remark:  In  the  above  proof,  we  can  see  that  if  \\T  — Th\\  < jh  then  ||Q  — Qh\\  < C-)h- 
This  gives  us  the  convergence  rate  of  the  eigenfunctions. 

The  next  theorem  is  proved  in  [27]. 

Theorem  4.4.4  Let  Aq  be  an  isolated  eigenvalue  of  T with  eigenfunction  uq.  Let 
{A/i}  be  a sequence  of  eigenvalues  of  Th  such  that  \h  —>■  Aq  and  {u/i}  C 14  the 
corresponding  eigenfunctions  that  converges  to  Uq.  Then  the  rates  of  the  convergence 

|A,-Ao|  < Ct{h)\ 


are 


< Ce{h), 


(4.31) 

(4.32) 


57 


for  sufficiently  small  h,  where  the  gap  e(h)  depends  on  the  choice  of  the  approximate 
subspaces  Vh- 

4.5  The  Subspaces  Vh 

Let  / :=  [0,a]  where  a is  the  core  radius.  We  define  the  approximate  spaces 
Vfi  as  follows.  Let  "P  be  a partition  of  the  closed  interval  I :=  [0,a]  defined  by 
0 = Xo  < xi  < . . . < xn  = a.  Set 

hi  :=  Xi  — x._i, 

h :=  max  hi, 
i 

min,  hi 
^ h ■ 

Let  Sh  be  the  finite  dimensional  space  of  continuous  functions  which  are  polynomials 
of  degree  1 on  each  subinterval  f :=  (x,_i,  Xi),  i — 1,2, . . . , N.  Define 

14:={u€  V:u|(0,a)e5,}.  (4.33) 

Thus,  Vh  consists  of  functions  in  V outside  the  interval  (0,a).  So,  Vh  is  infinite 
dimensional,  thus  the  essential  spectrum  of  Th  might  not  be  empty. 

We  now  show  that  the  essential  spectrums  of  T and  of  Th  are  identical.  This 
can  be  accomplished  if  we  can  show  that  the  difference  T — Th  is  a compact  operator 
(20] . 

Lemma  4.5.1  T — Th  is  compact. 

Proof:  We  want  to  show  that  for  any  bounded  sequence  in  H,  there  exists  a subse- 
quence whose  image  under  T — Th  is  a convergent  sequence.  We  consider  T — Th  on 
two  subintervals  (0,a)  and  (a,oo). 

1.  Let /„€  i/ with  ||/„||  = 1,  n = l,2, ... 
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We  want  to  show  that  the  sequence  {(T  — Th)fn}  is  precompact  in  H,  i.e.,  there 
exists  a subsequence  that  converges  (strongly)  in  H.  Denote 

Un  :=  Tfn 

y-hn  * — T^hfni  ^ — 1,2,... 

On  the  unbounded  subinterval  (a,oo),  u„  and  are  both  in  the  same  space  V, 
and  they  both  solve  the  non-homogeneous  equation  (4.17).  Hence,  the  difference 
Wn  :=  Un  — Uh„  satisfies  the  modified  Bessel  equation 
1 

— {rw'J'  + (a  - k‘^nl,)wn  + = 0,  (a,  oo). 

r r 

Since  this  equation  is  of  the  limit-point  at  infinity,  there  is  only  one  fundamental 
solution  which  is  denoted  by 

Ua{r)  = Km  . 

Therefore,  ru„(r)  must  be  a constant  multiple  of  Ua(r),  i.e., 

"Wnir)  = CnUa{r),  r > a.  (4.34) 

2.  By  the  proposition  4.2.2,  for  every  n, 

||^n||//‘(0,oo,r)  — ((^  T ^n)  = (/nj  rf„) 

< ll/n|lll«ni|  (II  • II  norm  in //■) 

= ll/n||||7^/n|| 

< ll/nfmi  = lirii  (||/„||  = 1). 

Therefore,  on  the  subinterval  (0,a),  we  obviously  have  ||^^n||^i(o  a r)  — ll^ll 
n.  Similarly,  ||w/,„||Hi(o,a,r)  ^ Hence, 


||^n||//l(0,a,r)  — H^n  ||i/*(0,a,r)  ^ 2||T||. 


(4.35) 
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3.  Here  we  show  that  the  constants  c„  in  (4.34)  are  bounded.  We  have 

{A  + a)(u„  - UhJ  (u„  - Uhjrdr  = 0.  (4.36) 

Integrating  (4.36)  by  parts  and  since  Ua{r)  > 0,  u'^{r)  < 0,  we  have 

(a  — k^rili)w^rdr  + I —dr  = —ac^u'^{a)ua{a) 

= acl\u'^{a)ua{a)\.  (4.37) 


However,  from  (4.35),  we  get 


cl< 


2\\T\\ 


(4.38) 


«K(«)wa(a)r 

4.  Finally,  by  (4.35),  iu„  is  a bounded  sequence  in  i/^(0,a,r).  For  bounded 
interval  (0,a),  H^{0,a,r)  is  compactly  embedded  in  F^(0,  a,r).  Hence,  there  exists 
a subsequence  of  labeled  again  by  {iy„},  that  converges  in  i^(0,a,r)  to,  say 

w 6 L^(0,a,r).  For  this  sequence,  we  also  have 


sup  |tn„(r)|  = sup  |c„||ua(r)|  < C\ua{r)\  Vr  G (a,  oo). 

n n 

Hence,  we  pick  a subsequence  of  {c„}  such  that  c„  — > Cq,  so 

|iin(r)  - u/»„(r)|  = |tu„(r)|  A |u;(r)|  in  1^(0, a, r),  (4.39) 

\un{r)  - Uh„{r)\  = k„(r)|  A co|ua(r)|  Vre(a,oo).  (4.40) 


These  imply  that  {io„}  contains  a subsequence  that  converges  in  H,  i.e.,  {T  — Th)un 
consists  of  a convergent  subsequence  in  H.  Thus,  T — T^,  is  compact  as  desired.  ■ 
Now  we  show  that  HT  — T)i||  < Ch^  for  some  s > 0 for  Vh  chosen  above.  On 
the  subintervals  /,-,  i = 1, 2, . . . , N ^ the  standard  approximation  theory  applies.  For 
the  interval  /q  = (0,  xj)  that  contains  the  singularity  xq  = 0,  we  need  a Poincare-type 
inequality. 
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Lemma  4.5.2  Let  u € H^{0,h,r)  with  u{h)  = 0.  Then 


f u^rdr  < h f (u')‘ 
Jo  Jo 


rdr. 


(4.41) 


Proof:  As  usual,  we  first  show  that  the  inequality  holds  for  functions  in  C^(0, /i]. 
Then  by  completion,  it  holds  for  functions  in  h,r).  If  u{h)  = 0,  we  see  that 

\y/ru{r)\  = \y/r  j u{t)'dt\ 

rh 

< I / \/tu'{t)dt\ 

rh 

< / Vi\u'{t)\dt  = 

Jo 

This  and  Cauchy-Schwartz  inequality  imply  that 

[ \u\hdr  < h||\/^w||ioo(o,;,) 

JO 

< h\\y/ru'\\li^o,h) 

Hence,  the  Poincare-type  inequality  holds.  ■ 

Lemma  4.5.3  For  each  u , there  exist  a constant  C = C{h,p)  and  w E Vh  such 
that 

||u'  - u;'||L2(0,a,r)  < Ch\\u"\\L2^o,a,r)  (4.42) 

where  h,  p are  defined  in  the  beginning  of  the  section. 

Proof:  Let  w{r)  be  a linear  polynomial  on  Iq  with  iy(a;i)  = u(xi)  and  iy'(a;i)  = u'{xi). 
Applying  lemma  4.5.2  to  the  difference  u'{r)  — w'{r),  we  get 

f>h  rh 


pti 

I |u'  — w'l^rdr  <h^  |u'^|rdr. 

Jo  Jo 
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On  each  i > 1,  let  w{r)  interpolate  u{r)  at  end  points  a;,}.  This  type  of 

interpolant  is  known  to  satisfy  [12] 


\u'  - tn'||L2(/i)  < 


(4.43) 


So, 


^^rdr 


' - r («'  - 

Jxi^l 

< Xi  j {u'  — w'Y 

Jxi-\ 

< Xih\  f [u")^dr 

Jxi-i 

< Xih^ f [u")^rdr 

J xi-i 

< C'(/i,p)/i^||u"||/,2(/,^r),  ^ = 2,3,. 


Now,  summing  over  j,  we  get 


^1|L2(0,a,r)  < Ch\\u"\\u2(Q^a,r)- 


(4.44) 


This  completes  the  proof.  ■ 

Now  we  are  ready  to  prove  the  error  estimate  for  the  non- homogeneous  problem. 
Recall  that  the  symmetric  V-elliptic  and  continuous  bilinear  form  b{u,v)  is  defined 
as 


f°°  /’°°  fm^  1 

6(u, u)  :=  a(u, u)  + a(u, u)  :=  J u'v'rdr J ^ — k^n^{r) -{■  a>  uvrdr. 

Theorem  4.5.4  Let  u E V and  Uh  € 14  be  the  solutions  of  the  non-homogeneous 
equations 


b{u,v)  - (f,v)  WvGV, 

b(uh,Vh)  = if,Vh)  Vuft  G 14, 
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where  f E H 


||w  — 

< C{Kp)h\\u"\\H, 

(4.45) 

||t^  - uh\\n 

< C{Kp)h^\\u"\\H. 

(4.46) 

Proof:  By  the  well-known  Cea’s  lemma  [6],  we  obtain 

b{u  -Uh,u  - Uh) 

= C inf  b(u  — Vh,u  — Vh) 

= C inf  y — Vh'Y  + - k'^n'^ir)  + {u  — rdr.  (4.47) 

There  are  two  cases: 

(i)  m = 0.  In  this  case,  the  coefficients  of  b(u,v)  are  bounded,  i.e.,  0 < 

Pn2(r)  < < a. 

(ii)  m / 0.  In  this  case,  we  have  seen  in  the  proposition  (4.2.1)  that 


r°°  ■m? 

I ~y(^  ~ Vh)^rdr  < 2 • / (u'  — v^Yrdr 

Jo  ^ Jo 


< oo. 


Hence,  by  the  definition  of  the  subspace  14  (14  = y on  (u,oo)),  the  remark  above, 
and  lemma  4.5.3,  there  exists  rn  € 14  such  that  the  infimumon  the  right  side  of  (4.47) 
can  be  reduced  to 

b{u  — Uh,u  — Uh)  < C inf  f {{1  + 2m^){u' — Vh)^  + {a  — k^n^(r)){u  — Vh)^}  rdr 

Vhev^  Jo 

< C f {{I  + 2m^){u  — w')^ {a  — k^n^{r))(u  — rdr 
Jo 

Therefore, 


b{u  -Uh,u-  Uh)  < C{h,p,m'^,a)  h?  ||w"||i2(o,a,r)- 
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The  continuity  of  the  bilinear  form  also  yields 

< C • h(u  - Uh,u  - Uh)^ 

< C'(/i,/?,m^a)  /i  ||n"||£,2(o,a,r)- 


(4.48) 


This  proves  the  first  part  of  the  theorem.  For  the  i/-norm,  we  use  Aubin-Nitsche’s 
duality  argument  [6].  Since  V C H is  dense  in  H,  we  have 

||«  - Uh\\H  < M\\u  - u/jIIv  sup  inf 

geH'l'heVh  ll^rllH 

where  for  any  g E H , x/^g  £ V \s  the  unique  solution  of  b{v,rpg)  = (g,v)H.  Using  the 
similar  argument  as  above,  we  can  easily  find  that 


inf  llV’a  - ^h\\v  < inf  b{xpg  - xl^h,  - i>h) 

i>h^Vh 

— w'Yrdr  + f {a  — k^n^){xpg  — w)^rdr 

Jo 

< C{h,p,m\a)  {h^||V’5llL(0,a.r)  + Ili2(0,a,r)} 

< C{h,p,m'^,a)  hJ  HV’^  |li2(o,a,r) 

< C{h,p,wJ,a)  \\g\\H  by  lemma  4.3.1. 

Thus,  combining  this  inequality  with  (4.48),  we  finally  have 

||n  - UhWn  < M h'^  ||w'1|L2(0,a,r)- 

This  completes  the  proof.  ■ 

Since  u = T f and  Uh  = Th.f,  we  just  show  that  with  this  particular  choice  of  14, 

\\{T-T,)f\\H  = 0{h^),  \\{T-n)f\\v  = 0{h)  (4.49) 

hence,  IIT-r,, II  = 0{h).  (4.50) 

So,  by  the  theorem  4.4.4,  the  rate  of  convergence  of  the  eigenfunctions  is  0{h)  and 
that  of  the  eigenvalues  is  0{h'^)  since  the  gap  c{h)  < ||T  — T/i||  by  definition. 
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4.6  Numerical  Experiments 

We  use  the  construct  in  the  previous  section  to  compute  the  propagating 
fields  in  circular  fiber.  We  now  consider  the  discretization  of  (4.8).  Since  we  know 
the  fundamental  solution  outside  of  the  interval  (0,a),  by  integration  by  parts,  the 
solution  Uh.  of  the  variational  formulation  satisfies 

b{uk,  Vh)  - \{uh,  Vh)  = 0 (4.51) 


where 


b(uh,Vh)  :=  / u'f^Vh'xdx+  / (a  - H -)uhVhxdx  - acau'^{a)vh{a),  (4.52) 

Jo  Jo  ^ 


and  Ua{x)  = Km{Wxfa),  and  A = a — j3‘^. 

The  finite  dimensional  Sh  has  a basis  {V’j}j=i,2,...,iv  where 
{x  - Xi) 


= (X.-.0 

for  1 < i < — 1, 


, X e (0,Xi),  if  m = 0 (see  the  next  remark  for  m > 1), 


4)j{x) 


and 


iPn{x)  = 


{x  Xj-i) 

{xj  ^j-l) 
{x:-Xj+i) 
(Xj-Xj+i) 

{x  - xn) 


X G (^j— i5^j)j 
X e {xj,xj+i), 

X e (xn-1,Xn). 


(x;v  — a^Af-l) 

Remark:  If  m > 1,  the  function  i/>o  is  excluded  from  the  basis  of  the  approximating 
subspaces  since  u(0)  = 0. 


Then,  each  G Sh  can  be  expressed  as 

N 

Uh  = Y^  UjXpj 
j=l 


where  uj  = Uh{xj). 


(4.53) 
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Denote  the  matrix  A = (aij)  where 

a,j  = J tpl{x)il^j{x)xdx  + j (a  — k^n^{x))ipi{x)ipj{x)xdx 
and  the  matrix  Bq  = (bij)  where 

6,  = /*(x)^,(x)x<(x. 

Denote  U = {uj).  Then  the  matrix  formulation  of  (4.51)  is 


(4.54) 


(4.55) 


AU  - XBoU  + (0, . . . , aCauMf  = 0.  (4.56) 

Eigenfunctions  Uh  and  their  derivatives  u),  are  continuous  on  (0,  oo),  so  un  •= 
Uhi^N)  = Uh{a)  = CaUa(a)  — CaKm(W).  Hence,  the  constant  Ca  is 


Un 

Ca  — Z T • 
Ua{a) 

Substituting  this  in  (4.55),  we  get 


{A-B{X))U  = 0 
where  B(A)  :=  AJ5o  — diag(0, . . . , /(A)) 

Ua(a)  Km{W) 

with  W — ayj a — X — 


Therefore,  our  task  is  to  solve  for  the  eigenvector  U and  the  eigenvalue  A.  For  this, 
we  apply  the  inversion  algorithm  with  shift  strategy  as  in  [17]. 

Algorithm;  Let  A and  5(A)  be  Hermitian  in  a neighborhood  of  Aq  where  Aq  is  an 
isolated  eigenvalue  of  {A  — B{X))U  = 0.  Assume  that  the  derivative  B'{X)  is  positive 
definite  in  Given  Xh  € Ue  and  Uh  with  Xh  / Aq  and  B'{Xo)Uo  / 0.  Define 


{A-B{Xh))V  ;=  B'{Xh)Uh, 


(4.57) 
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'/i+i 


+ 


V^B'{\H)Uh 

V^B'(Xf,)V 


Uh+,  := 


V 

iiv'ir 


Then  for  sufficiently  small  e,  Xh  ->  Aq  of  order  C?(/i^)  and  Uh  Uo- 


(4.58) 

(4.59) 


Example:  Since  exact  solutions  for  the  circular  fiber  with  truncated  parabolic 
index  profile  are  possible,  we  can  check  the  accuracy  of  our  finite  element  solutions. 
We  recall  that  the  truncated  parabolic  profile  n[r)  is  defined  by 


n^(r) 


n: 


(l-2Ag) 


0 < r < a, 
r > a, 


where 


Ur 


= 1.48//m,  rici  = 1.46/xm,  a = 31.25//m, 


A = 0.85/um(wavelength),  A :=  - ^1 ^ ) • 


In  chapter  3,  we  have  seen  that  the  exact  solutions  for  TE  case  are  in  terms  of  the 
Whittaker  functions.  The  propagation  constants  are  found  by  solving  the  eigenvalue 
equation 

2V^~~-  - 1 = (4.60) 


M«,,(l/)  ^ ■■  Km{W) 

where  Mk,^(z)  is  the  Whittaker  function  and 

V = akyjnl^  - n%  U = ayjk^nl^  - /?2,  W = ayj - k'^nh, 


m 


K = 


AV'  ^ 2 


The  first  two  largest  propagation  constants  /?  for  m = 0, 1 are  computed  by  the  finite 
element  method.  They  are  compared  to  the  exact  eigenvalues  found  in  (4.60).  The 
corresponding  modes  are  well  confined  inside  the  core  of  radius  a.  The  following 
tables  are  self-explanatory. 
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Figure  4.1;  First  mode  for  m = 0,/i  = a/128 


Figure  4.2;  Second  mode  for  m = 0, /i  = a/512 


68 


Figure  4.3:  First  mode  for  m = l,/i  = a/512 


Figure  4.4:  Second  mode  for  m = l,/i  = a/512 
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Table  4.1:  Approximations  of  the  largest  /3  for  m = 0. 


N 

Finite  element 

Error 

Rate 

Exact 

8 

10.9347148 

1.7e-4 

16 

10.9348452 

4.5e-5 

1.323 

32 

10.9348785 

1.2e-5 

1.907 

64 

10.9348872 

2.9e-6 

2.048 

128 

10.9348894 

7.0e-7 

1.421 

256 

10.9348900 

l.Oe-7 

1.946 

512 

10.9348901 

2.5e-8 

2.000 

10.9348901 

Table  4.2:  Approximations  of  the  second  /?  for  m = 0. 


N 

Finite  element 

Error 

Rate 

Exact 

8 

10.9228920 

1.5e-3 

16 

10.9239979 

4.0e-4 

1.901 

32 

10.9242895 

l.Oe-4 

2.000 

64 

10.9243675 

2.7e-5 

1.889 

128 

10.9243874 

6.8e-6 

1.989 

256 

10.9243925 

1.7e-6 

2.000 

512 

10.9243937 

5.0e-7 

1.765 

10.9243942 

Table  4.3:  Approximations  of  the  largest  /?  for  m = 1. 


N 

Finite  element 

Error 

Rate 

Exact 

8 

10.9293260 

3.2e-4 

16 

10.9298407 

1.9e-4 

0.752 

32 

10.9297378 

9.4e-5 

1.015 

64 

10.9296700 

2.6e-5 

1.854 

128 

10.9296502 

6.4e-6 

2.022 

256 

10.9296451 

1.3e-6 

2.299 

512 

10.9296438 

3.0e-7 

2.115 

10.9296438 
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Table  4.4;  Approximations  of  the  second  /3  for  m = 1. 


N 

Finite  element 

Error 

Rate 

Exact 

8 

10.9164070 

2.7e-3 

16 

10.9187833 

3.6e-4 

2.906 

32 

10.9191934 

5.1e-5 

2.819 

64 

10.9191672 

2.5e-5 

1.028 

128 

10.9191491 

6.7e-6 

1.899 

256 

10.9191441 

1.7e-6 

1.979 

512 

10.9191428 

4.0e-7 

2.086 

10.9191424 

Table  4.5:  Error  in  weighted  L^-norm,  m = 0. 


N 

Error 

Rate 

128 

.0172 

256 

.0073 

1.231 

512 

.0030 

1.274 

1024 

.0012 

1.304 

2048 

.0004 

1.328 

4096 

.0001 

1.346 

Table  4.6:  Error  in  L^-norm,  m = 0. 


N 

Error 

Rate 

128 

.0096 

256 

.0042 

1.201 

512 

.0017 

1.251 

1024 

.0007 

1.286 

2048 

.0002 

1.313 

4096 

.0001 

1.334 

CHAPTER  5 

NON-CIRCULAR  FIBERS  AND  DTN  MAPS 


In  this  chapter,  we  consider  weakly  guiding  fibers  with  arbitrary  core,  i.e., 
not  necessarily  circular.  Hence,  we  need  to  solve  Helmholtz  eigenvalue  equations 
governing  the  guided  modes  for  the  whole  plane  IR^.  Combining  the  standard  finite 
element  method  and  Dirichlet-to-Neumann  map  (DtN),  we  can  reduce  the  unbounded 
computational  domain  to  a bounded  one,  in  particular,  to  a circular  disk  of  radius 
R.  We  use  this  method  to  the  study  the  propagating  evanescent  wave  energy  that 
exists  outside  the  core  of  the  square  and  circular  fibers. 


5.1  Dirichlet-to-Neumann  Map 


Since  the  geometry  of  the  core  is  not  circular,  polar  coordinate  transformation 
is  no  longer  helpful  as  in  Chapter  3 and  Chapter  4.  We  now  consider  the  2D  scalar 
Helmholtz  equations  for  the  longitudinal  components  and  of  the  electric  and 
magnetic  fields  of  propagating  waves  in  the  fiber: 


(Po) 


and 


AE,  -k  (Pn2  - p^)E,  = 0 in  IR^ 

Ez,  are  continuous, 

on 

E,  e 

AH,  + - f3)H,  = 0 in  IR^ 


Rzt  — are  continuous, 
2 


n“ 


dn 


H,  £ 


(5.1) 


(5.2) 
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Vu  • Vu  dxdy  + [ -^uv  dxdy  — f uv  dxdy  (5.4) 

Jb?  Jb? 


To  ease  the  notation,  we  denote  u = or  The  variational  formulation 
of  (5.1)  is 

a[u,v)  :=  — / Vu  ■ Vu  dxdy  + / k^n^uv  dxdy  — (3^  I uv  dxdy  (5.3) 
Jb^  Jb^  Jb'^ 

and  that  of  (5.2)  is 
a{u,v):=  I ^ 

Jb^  ^ 

Remark:  In  the  variational  formulation  of  H^,  k^  plays  the  role  of  eigenvalue. 

We  concentrate  on  the  problem  of  Ez,  the  problem  can  be  treated  identi- 
cally. If  we  discretize  (5.3),  we  face  the  difficulty  of  the  unbounded  domain.  One  of 
the  widely-used  techniques  is  to  use  compact-support  finite  elements  to  approximate 
the  unbounded  discrete  problem.  This  is  known  to  produce  many  spurious  solutions. 
We  hence  adapt  a method  known  as  Dirichlet-to-Neumann  map  used  in  exterior  and 
scattering  problems.  We  show  later  that  this  method  eliminate  the  undesired  so- 
lutions. There  are  three  main  steps  in  using  DtN  methods.  First,  an  appropriate 
artificial  boundary  F/i  needs  to  be  introduced  so  that  its  interior  contains  the  core 
region  Cl.  Second,  a boundary  condition  cissociated  with  F/j  has  to  be  determined  in 
such  a way  that  a Dirichlet  condition  is  related  to  a Neumann  condition.  Third,  a 
new  reduced  ’’interior”  problem  needs  to  be  solved. 

(i)  Artificial  boundary  Fh;  A natural  choice  for  Fh  is  a circle  of  radius  R where  R 
is  sufficiently  large  so  that  Cl  is  contained  in  the  interior  of  F^.  The  reason  of  this 
choice  will  be  clear  later. 

(n)  The  DtN  map:  Let’s  denote 

Cli  ;=  Bft  the  disk  of  radius  /2-the  inferior  domain 


and 


Cle  = IR^  \ Br  the  exterior  domain. 
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Hence,  the  Helmholtz  problem  in  (5.1)  is  equivalent  to  the  following  eigenvalue  prob- 
lem: 


Let  Ui  :=  and  Ug  = w|n^.  Then 

Aui  -t-  {k^n^  — P^)ui  = 0 in  Qi, 
Aue  + {k'^n^  — I3‘^)ue  = 0 in  fie, 

Ui  = Ug  on  Tr, 
dui  dug 

on  1 R. 

k C/Tl  C/Tl 


(5.5) 


Assume  that  we  can  solve  for  Ug  analytically.  We  want  to  construct  a map  that 
relates  Ui|r  and  ^|r-  The  next  theorem  is  a standard  theorem  for  trace  mappings 


[3] 


Theorem  5.1.1  Ifdfl  is  a smooth  boundary  offl  then  there  are  continuous  (bounded) 
linear  mappings  (called  traces) 

T -.u^ulan:  (5.6) 


For  our  problem,  the  Hilbert  space  is  V :=  H^{Br)  with  the  smooth  boundary  Fh- 
Hence,  € H2(TrY  We  then  define  the  map 

Tr  : //'/'(Fh)  ^ /7-'/"(Fh) 


by 


TH(u,|rJ  = ^ 


Consequently,  the  boundary  condition 


dui 

dn 


dug 

dn 


Fr 


becomes 


TR{ui\r,^) 


dug 


dn 


Fr 


(5.7) 


(5.8) 
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Before  we  construct  Tr  explicitly,  we  consider  the  following  auxiliary  problem  in  the 

exterior  domain  fig: 

Given  k'^nh  < and  g 6 

find  Ue  G H^{Vle)  such  that 

f Aug  + - /?^)ug  = 0 in  fig, 

\ Ue  = g on  Tr. 

Lemma  5.1.2  The  auxiliary  problem  (5.9)  has  a unique  solution  Ug  G H^{ile). 


Proof:  Let’s  set  u;  = itg  — ^ where 


9 - 


g on  Tr 
0 elsewhere. 


Hence,  u;  = 0 on  Ph-  Hence,  we  obtain 


Ao)  + {kin'll  — = 0 in  fig, 

a;  = 0 on  Ph- 


(5.10) 


Let  H = and  V = {to  e : a;  = 0 on  P/j}.  Then  V C H is  a Hilbert 

space.  The  variational  formulation  of  (5.10)  is 

a(iu,u):=  / Ww  ■ Vvdxdy {j3^  — k^nh)  I wvdxdy. 

Jug  Jne 

Since  (3^  — k^n(  > 0,  a(-,  •)  is  coercive  and  bounded.  By  Lax-Milgram  theorem,  (5.10) 
has  unique  solution.  This  proves  the  lemma.  ■ 

With  this  existence  and  uniqueness  result,  we  can  find  the  solution  Wg  in  (5.9)  explic- 
itly. In  polar  coordinates,  the  exterior  problem  becomes 

d Ug  1 dug  1 d Ug  . 2 I 2 2 \ n • o 

, + + 

. Ue\rn  = g{d)  on  Pfl. 

Using  the  Fourier  series  representation  of  Ug,i.e., 


u, 


m=0 


(5.11) 
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we  obtain 


We  denote 


dr‘2  i>  dr 


n 


cl) 


u 


= 0. 


a = 


(5.12) 


For  each  n = 1, 2,  • • •, 

= AnKn{ar)  + Bnln(ctr). 
However,  only  Kn{ar)  is  an  function.  Thus 

CX> 

= y^,  AnKn{ar)e‘’^^ . 


(6.13) 


(5.14) 


n=0 


On  the  boundary  F/?,  we  also  express  the  boundary  function  g{6)  in  Fourier  series, 


i.e.. 


gW  = 


n=0 


Then,  the  condition  Ug  = g on  T r yields 


A„A'„(afl)  = jW  n = 0,1,2, 


(5.15) 


or 


An  — 


il”) 


KniaR) 


Therefore,  the  solution  Ug  is  of  the  form 


We  note  that 


dug 


dn 


°° 

= E rT-pi  ■ A'»(a>-)^"'’. 

dug 


Tr  dr 


n Jn)Jn6 

— / ^ ITT riT  9 ^ 


rfi  ^ Rn{aR) 


(5.16) 


(5.17) 


(5.18) 
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Thus,  we  can  define  the  map  Th  : H^{Tr)  H ^(r;^)  as  follows 


E“ 


An)Jne 

Kn{aR) 


oo 

n=0 

In  the  next  lemma,  we  study  some  of  the  properties  of  Tr. 


Lemma  5.1.3  Tr  : H‘i{TR)  »->  H ^(r;?)  defined  by 

rp  ± J^(n)  ine 

is  a continuous,  symmetric,  and  linear  operator. 

Proof:  The  linearity  of  Tr  is  obvious.  We  note  that  if  ^ G H2(Tr), 

< Tr</>,  ^ > :=  TR(j)-^drR 
Jvr 


/•27r  / °°  o°  \ 

= ^ Rde 

\n=0  n=0  / 


= 27Ti?  ^ 


n=0 


n=0  ' 


By  the  definition  of  dual  norm,  we  have 


W'^R^WR-hr  ) ■=  sup  |<Tr0,^> 
" ll^ll  i =1 

H2I 


Urfi) 


< sup 
ll'I'll  1 =1 


2ttR  Q 


H2 


(r-fi) 


n=0 


K'(c^R) 


(5.19) 
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Since  K'^{aR)  < 0 and  Kn{aR)  > 0,  we  can  write 


K{cR) 


Kn{aR) 

Using  the  following  identity  in  [1] 


-Kic^R) 

KniaR) 


zK'„{z)  = -nKn{z)  - zKn-i{z), 


we  get 


-Kjc^R) 

Kn{aR) 


_ , Kn-i{aR) 

aR  Kn{oiR) 

n 


Thus,  for  sufficiently  large  n,  we  have 

-^  + 1 < (n^  + 1)'/^ 

aR 

In  fact,  the  inequality  holds  for  n > 2aRI{a^ R?  — 1).  Thus,  for  n 1,  we  obtain 


E 


n=0 


K{aR) 


|^W|2  < . |^(n)|2_ 


n=0 


KniaR) 

Hence,  by  Cauchy-Schwarz’s  inequality,  we  obtain 

K(aR) 


E 

n=0 


K4aR) 


So,  Tr  is  bounded  in  The  symmetry  of  Tr  is  clear  from  definition.  ■ 

Now,  we  are  ready  to  define  the  interior  problem. 
iiii)  The  interior  problem  in  Q, 

Find  i/3\ui)  eR*x  H\Br)\{0} 

Au{  + (k^n^  — = 0 in  Br^ 

{Pi)  { (5.20) 


dn 


= 7H(w,|rR)  on  Tr. 


We  need  show  that  the  new  problem  is  equivalent  to  the  original  problem  (Po)- 
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Lemma  5.1.4  (/3^,u)  € /?*  x \ {0}  is  a solution  of  (Po)  iff  (P^,Ui)  € R*  x 

\ {0}  is  a solution  of  (Pi). 


Proof:  Suppose  (f3'^,u)  is  a solution  of  (Po)-  Let  u,-  :=  w|b^-  Then  Ui  ^ 0;  otherwise, 


by  the  continuity  conditions  Ug  = Ui  and 
due 


dug  dui 


on  Pk,  we  would  get  Uelrn  = 0 


and 


dn 
see  that 


dn  dn 

0.  It  would  imply  that  Ue  = 0,  so  u = 0.  By  the  definition  of  Tr,  we 


TR(u,|rH)  = Tfi(nelrH) 

dug 


dn 

dui 

dn 


r« 


Thus  u|n,  also  satisfies  the  boundary  condition  in  (P,).  Consequently,  let  (/?^,u,)  be 
a solution  of  (P,  ).  Let  Ug  be  the  unique  solution  of  the  exterior  domain  with  the  fixed 
/3^  with  boundary  condition  Ue|rH  = Hence,  we  see  that 

dug 


dui 

dn 


rfi 


= TH(u.lrH)  = Pfl(uelrH) 


dn 


Tr 


Thus,  Uelrn  = ■“ilrn 


dup 


dn 


Tr 


dui 

dn 


are  satisfied.  The  function  u whose 


restrictions  on  fl,-  and  fig  are  U{  and  Ug  is  a solution  of  (Po).  This  u is  unique  because 
Ug  is  unique  (by  Lax-Milgram  for  the  exterior  problem).  ■ 

5.2  Variational  Formulation  of  Interior  Problem 


We  recall  that  H = L^(Br)  is  a Hilbert  space  with  the  usual  inner  product 
and  L^-norm  and  V :=  H^{Br)  C H he  a subspace.  Let  Ua  be  an  arbitrary  constant 
such  that  Ua  > Ugo-  The  bilinear  form  6(-,  •)  defined  on  V x V is 


b{u,  v) 


Vu  • Vu  dxdy  + 


— n^)uv  dxdy  — 
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= / Vu  • Vu  dxdy  + k^  (n^  — n^)uv  dxdy 

J Bn  J Bn 

^ Kn{aR) 


n=0 


(5.21) 


Note  that  we  have  introduced  the  term  k^n'^[u,v)  into  the  variational  form. 


Lemma  5.2.1  The  bilinear  form  b(-,-)  is  symmetric,  bounded,  and  coercive,  i.e., 
there  exists  Ci , C2  > 0 such  that 


|6(u,u)|  < Cillullv  ||u||v, 
b{u,u)  > C2\\u\\v- 

Proof:  The  symmetry  of  b{u,v)  follows  from  its  definition,  i.e., 


b{u,  v)  = b{v,  u). 


Since  the  linear  map  Tr  is  bounded  in  H 2(Fh),  we  have 


Hence,  by  the  Schwarz’s  inequality,  we  have 


(f  TR{u\rn)vdrR 
Jrn 


< C'||u||^i/2(rH)||i’||//i/2(rH) 

< C'||w||//‘(Bh)II^IIh‘(Bh)- 

The  last  inequality  follows  from  the  trace  theorem  mentioned  earlier.  Thus,  we  have 

|6(u,u)|  < ||Vu||//||Vu||// + - n^i)||u|li/||u||// + CIIuIIkIIuIIi^ 


< C||u||v||u||y. 
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Hence,  6(u,  v)  is  bounded  in  V x K.  To  show  that  b(-,  •),  we  notice  that  — > 0, 
we  get 

b(u,  u)  > ||Vu||^  + P(nl  - n^)\\u\\l  > C\\n\\l. 

Thus,  6(-,  •)  > 0 for  non-zero  u G V.  ■ 

Lemma  5.2.2  For  each  f G L‘^{Br),  there  is  a unique  solution  u E V to  the  non- 
homogeneous  variational  equation 

b{u,v)  = {f,v)  VuGK  (5.22) 


Moreover,  u G H^{Br)  and 


\u\\H^BR)  < C\\f\\L2(BR)- 


(5.23) 


Proof:  The  existence  and  uniqueness  of  the  solution  u G V is  the  standard  con- 
sequence of  Lax-Milgram  theorem.  The  regularity  of  u is  also  a standard  result  in 
potential  theory  [22]  where  the  non-homogeneous  problem  is  viewed  equivalently  as 


Au,  + k‘^{nl  - n^)u,  = f 
Aue  + k^{nl  - nh)ue  = 0 


Ui  = Ue 

dui  due 


dn  dn 

dufi 


in  Br, 
in  IR^  \ Br, 
on  Fr, 

on  Tr, 


a/t  — > 0,  a/t  Ue  — > 0 as  r ^ oo. 

or 

This  completes  the  proof.  ■ 

By  a representation  theorem  for  bounded  symmetric  bilinear  forms,  there 
exists  a unique  bounded  self-adjoint  operator  B defined  on  D(B)  C V such  that 


b{u,v)  = (Bu,v)  Vu,uGK 
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and  B has  the  same  bounds  as  6(-,  •).  Since  0 is  a point  in  the  resolvent  set  of  5,  the 
resolvent  T :=  is  well-defined,  i.e.,  T is  bounded. 

Lemma  5.2.3  T is  a compact  operator  in  H . 

Proof:  Since  H^{Br)  is  compactly  embedded  in  77  = L'^{Bji), 

T:H^  D{B)  C V = H\Br)  ^ H. 

Therefore,  T is  compact.  ■ 

This  implies  that 

Theorem  5.2.4  [20]  The  spectrum  of  B consists  of  entirely  of  isolated  eigenvalues 
with  finite  multiplicities  and  72^(77)  :=  {B—Q~^  is  compact  for  every  ^ in  the  resolvent 
set  of  B . 

5.3  Finite  Element  Approximation 

This  section  is  very  similar  to  section  4 in  Chapter  4.  Suppose  that  V/,  is  a 
family  of  finite  dimensional  subspaces  of  V that  satisfies  the  approximation  properties. 
Approximation  Properties  of  14:  For  any  u € H\Bp),  there  exists  an  interpolant 
7T/,u  € 14  such  that 

\W  - ^hu\\Hi(Bn)  < Ch‘~‘\\u\\H‘{BR),  S > /,  (5.24) 

where  h is  the  mesh  size.  Since  we  mainly  work  with  triangular  linear  elements,  h is 
the  maximum  of  the  diameters  of  triangles  in  a triangulation. 

As  in  Chapter  4,  we  first  study  the  discrete  inhomogeneous  problem: 
let  / G 77,  find  Uh  G 14  such  that 


6(u,,u'‘)  = (/,u'*)  Vu'^G  V4. 
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Let  {Bh,  D{Bh))  be  the  unique  self-adjoint  operator  associated  to  6(-,  •)  : V/i  x \4  C. 
Let  Th  :=  Bf^^.  Then  is  compact,  in  fact,  Th  is  of  finite  rank.  Then  a unique 
solution  Uh  to  the  inhomogeneous  problem  is  of  the  form 

UH  = T,/,  / e i/.  (5.25) 

We  have  the  same  convergence  rate  as  in  the  singular  Sturm-Liouville  case,  i.e., 
theorem  4.5.4. 


Theorem  5.3.1  Let  u and  Uh  be  the  unique  solutions  of  the  variational  problems: 


b{u,v)  = (/,u),  veV, 

(5.26) 

b{uh,v^)  = (/,u^),  v^eVh. 

(5.27) 

Then  they  satisfy 

llw-u/illv  < Ch\\u\\u2(^Bii), 

(5.28) 

\\u  — Uhlln  < Ch^\\u\\H2^Bii)- 

(5.29) 

Proof:  By  Cea’s  lemma,  we  have 

b{u  — Uh,u  — Uh)^  < C inf  b(u  — v^,u  — 2 

< C inf  llu  — 

< Cllii  - nhu\\v 

< Ch\\u\\H2^BR)- 

(5.30) 

But,  we  also  have 


b{u  - Uh,u  - Uh)'^  > C\\u-  Uh\\v, 
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thus, 

11^^  - Uh\\v  < Ch\\u\\ff2^BR)- 

Hence,  the  first  inequality  holds. 

Now  we  use  Nitsche’s  trick  as  follows.  Let  g ^ H . Let  Wg  ^ H he  the  unique  (dual) 
solution  of 

h{v,Wg)^{g,v),  veV. 

Then  we  have 

inf  h{wg  — rt)^,  Wg  — < Cb{wg  — nhW,  Wg  — nhw)^ 

w^eVh 

< C\\Wg  - TThWgWv 

< Ch\\Wg\\H^BR) 

< Ch\\g\\mBn), 

therefore,  \\wg  — w^\\v  < Ch\\g\\i,2^Bj^)  for  h sufficiently  small.  Thus,  by  Aubin-Nitsche 
Lemma,  we  get 

< C'||u-u;i||v  sup  inf  \\Wg-W^\\v 

< Ch'^\\u\\H^BR)- 

This  proves  the  second  inequality.  ■ 


Corollary  5.3.2  For  h > 0 sufficiently  small,  we  have 

m - T\\  < ch. 

This  corollary  implies  that  the  compact  discrete  operators  Th  converges  in  norm  to 
the  compact  operator  T.  Now  we  consider  the  variational  form  of  eigenvalue  problem: 

Find  (A,?i)  G DL*  x V such  that  b{u,v)  = \{u,v)  Vu  € V (5.31) 
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where  A = By  the  theorem  5.2.4  the  eigenvalues  A are  isolated  since  if 

fj,  is  an  eigenvalue  of  T then  A = A.  This  is  important  for  numerical  computations. 
The  associated  discrete  eigenvalue  problem  is 

Find  (Xh,Uh)  G IR*  X Vh  such  that  b{uh,v^)  = A/i(u/i,u^)  € 14.  (5.32) 

We  now  state  the  standard  theorem  in  finite  element  theory  on  the  convergence  rate 
of  the  eigenvalues  and  eigenvectors  for  the  operator  B with  compact  resolvent  T. 

Theorem  5.3.3  [12]  Let  Aq  be  an  isolated  eigenvalue  of  B with  a normalized  eigen- 
vector uq.  Then  there  exists  an  eigenvalue  Xh  of  Bh  with  a normalized  eigenvector  Uh 
such  that 

\Xh-X\<Ch'^ 


and 

||w;i  - u\\h  < Ch. 

5.4  Numerical  Experiments 

Now  we  present  some  numerical  experiments  of  the  finite  element  method. 
Let  us  recall  the  variational  form 

b{u,v)  = / VuWv  dxdy  / k^{n\  — n^)uv  dxdy  — (p  TnuvdYR 

JBfi  JBr  J^r 


where 


K'niar)  / . -Q 

Tau  = > a — — — L ^ 


n=0 


For  computation,  we  rewrite  the  series  expansion  of  the  boundary  operator  Tr  as 

Tru  = V Q cos(n0)  + sin(n6>)) 

^ Kn{aR) 
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u] 


u] 


where  u^^\r)  and  are  the  Fourier  coefficients  of  u defined  as 

‘fV)  ^ u{r,e)  d6, 

= — / u{r,6)cos(nd)  dO,  for  n > 1, 

^ Jo 

uf\r)  = 0, 

~ / u(r,  0)  sin(u0)  dO,  for  n > 1. 

Jo 

We  define  a triangulation  T in  the  bounded  computational  domain  Br  as  usual  and 
the  corresponding  partition  V of  the  circle  Tr.  Let  Vh  be  defined  as 

Vh  :=  K € V ; v%  € P\r),  reT,  e ^'(7),  7 € V}.  (5.33) 


We  consider  a step-index  circular  fiber  with  core  radius  a = 1,  core  refractive  index 
= 2,  and  cladding  index  rid  = 1-  Let  F^  be  a the  circle  of  radius  R = 2.  For 
the  circular  case,  there  exist  exact  analytical  solutions  in  terms  of  Bessel  functions 
Jn{Ur)  and  modified  Bessel  functions  A"„(VFr)  where 

U = W = ay^/?2  - k^nl 


and  the  associated  eigenvalue  equations  are  [5] 


fjML  + -^niU)  , KiW)  \ 

\UJn{U)  ^ WR\{W)J  \nl  UJn{U)  ^ WKr,{W)J 


4. 


(5.34) 


Therefore,  we  can  compare  our  results  with  the  exact  solutions.  As  we  can  see  from 
the  tables  that  with  only  coarse  meshes  we  can  obtain  good  approximations  to  the 
exact  eigenvalues.  The  errors  are  about  5 %. 
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Table  5.1:  FEM  eigenvalues  for  = 3.762625 


h 

FEM 

Error 

Exact 

.7789 

2.216394 

.4203 

2.253553 

0.037159 

.2908 

2.258267 

0.004714 

.1137 

2.258286 

0.000019 

2.2001 

Table  5.2:  FEM  eigenvalues  for  = 5.783529 


h 

FEM 

Error 

Exact 

.7789 

2.904700 

.4203 

2.953218 

0.048518 

.2908 

2.961140 

0.007922 

.1137 

2.965274 

0.004134 

.0916 

2.965319 

0.000045 

2.8930 

5.5  Square  Fiber  Versus  Circular  Fiber 

In  this  section,  we  apply  the  above  method  to  compute  the  evanescent  energy 
existing  outside  of  a fiber.  This  energy  is  important  in  optical  spectroscopy.  The 
cladding  of  the  fiber  is  removed.  The  core  is  then  used  as  sensor. 

We  consider  a square  fiber  of  side  s = y/2  and  a circular  fiber  of  radius  a = ^. 
The  circular  core  is  smaller  than  the  square  one.  The  evanescent  energy  of  each  mode 
excited  in  a fiber  is  the  T^-norm  of  the  solution  restricted  to  the  exterior  region.  Let 
rici  — 1 for  air  and  the  core  index  rico  = 1-5.  We  will  vary  the  wavelength  A and  observe 
the  change  in  the  evanescent  energy  in  each  fiber.  Since  the  lower  order  modes  are 
dominant,  we  concentrate  on  the  first  two  modes  {LPqi  and  LP02)  which  are  always 
excited  under  a normal  illumination.  In  the  figures,  we  plot  the  whole  intensity  of 
the  propagating  modes  over  the  computational  domain  and  also  the  portion  of  the 
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evanescent  intensity  outside  each  fiber.  We  use  the  meshes  with  size  h = .1128  for 
the  square  fiber  and  h = .1137  for  the  circular  one. We  observe  that  as  the  wavelength 
A increases,  more  evanescent  energy  exists  outside  of  each  fiber.  Though  the  square 
fiber  is  larger  in  core  area,  its  evanescent  energy  is  about  twice  as  that  of  the  circular 
fiber  for  both  considered  modes  {LPoi,  LPq2-  The  last  figure  shows  the  percentage  of 
evanescent  energy  as  a function  of  the  wavelength  A = 0.45/im, . . . , 2/im.  It  is  cleax 
that  the  square  fiber  releases  more  evanescent  energy  to  the  outside  region  than  the 
circular  fiber.  Hence,  it  seems  that  for  spectroscopy  application,  square  fibers  would 
be  a better  tool. 
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X 10  ^ 


X 10"'’ 


Figure  5.1:  LPq\  mode  in  circular  fiber,  A = 1.6/xm 
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X 10“^ 


Figure  5.2:  LPq2  mode  in  circular  fiber,  A = 1.6/im 
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Figure  5.3:  LPoi  mode  in  square  fiber,  A = 1.2/rm 
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Figure  5.4:  The  profile  of  the  evanescent  intensity  at  the  interface  of  LPoi  mode  in 
square  fiber 
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Figure  5.5:  LP02  mode  in  square  fiber,  A = 1.2/xm 
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Figure  5.6:  The  profile  of  the  evanescent  intensity  at  the  interface  of  LP02  in  square 
fiber,  A = 1.2/im 
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Figure  5.7:  The  evanescent  energy  vs  A.  Circular:o,  squareio.  Top:  LPoi  modes. 
Bottom:  L/02  modes. 
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